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FOREWORD 


This  document  contains  a  final  report  of  research  done  on  Contract  No. 
F49620-83-C-0096,  entitled  "Research  on  Turbine  Flowfield  Analysis  Methods". 

The  initial  effort  was  presented  in  an  interim  report  (Ref.  la),  and  a  re-.ent 
summary  was  given  in  Ref.  lb.  Technical  effort  was  curtailed  for  a  peri^  1  of  time 
when  the  Principal  Investigator  (Dr.  W.  J.  Rae)  became  a  member  of  the  culty 
in  the  Department  of  Mechanical  and  Aerospace  Engineering  at  the  Stat  livers ity 

of  New  York  at  Buffalo.  Arrangements  were  made  to  resume  the  effort  £  “le 

University  on  a  subcontract  basis  and  a  no-cost  time  extension  through  e  end 

of  1984  was  granted  to  facilitate  completion  of  the  program.  The  progr.  remained 
under  the  direction  of  Dr.  Paul  V.  Marrone,  Head,  Physical  Sciences  Depa.  ‘ment, 
Calspan  ATC. 

This  document  contains  the  completion  of  the  work  initiated  in  Ref.  1. 
The  results  presented  make  it  possible  to  generate  computational  grids  for 
cascades  of  high  solidity,  and  to  calculate  the  imcompressible  potential  flow 
on  such  a  grid. 


la.  Rae,  W.J.,  "Revised  Computer  Program  for  Evaluating  the  Ives  Transformation 
in  Turbomachinery  Cascades"  Calspan  Report  No.  7177-A-l  (AFOSR  Report  No. 
TR-83-1284)  July  1983 

lb.  Rae,  W.J.  (State  University  of  New  York  at  Buffalo)  and  Marrone,  P.V. 

"Research  on  Turbine  Flowfield  Analysis  Methods"  Calspan  Report  No.  7177-A-2 
April  1984 
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Sactlon  1 
IMTRODOCTIOM 

Some  of  the  most  Important  recent  advances  in  computational  fluid 
dynamics  have  been  made  possible  by  the  use  of  algorithms  which  solve  the 
equations  of  motion  In  boundary-conforming  coordinates  (see,  for  example,  Refs. 
2-4).  Thus,  the  development  of  methods  for  carrying  out  these  coordinate  trans¬ 
formations  has  received  considerable  attention  (for  example.  Ref.  5) . 

The  field  of  turbomachinery  flows  has  been  the  object  of  a  major  fraction  of 
this  attention,  partly  because  the  constraints  Imposed  by  Internal  flows  make 
the  application  of  certain  contemporary  algorithms  more  difficult  than  Is  the 
case  for  external  flows.  In  particular,  the  presence  of  non-orthogonallty  or 
shearing  In  the  computational  grid  tends  to  malce  the  algorithms  less  stable,  and 
the  solutions  are  more  sensitive  to  the  methods  used  for  treating,  In  the  com¬ 
putational  plane,  the  regions  corresponding  to  trailing  edges  and  the  points  at 
upstream  and  downstream  Infinity.  References  6^7  are  typical  of  one  ap¬ 
proach  to  the  resolution  of  the  latter  class  of  problems.  The  net  conclusion 
drawn  from  that  work  la  that  there  Is  a  need  for  Improved  grid -generation  tech¬ 
niques  for  turbomahclnery,  particularly  with  respect  to  the  treatment  of 
tralllng-edge  regions  and  to  the  achievement  of  orthogonal  grids. 

The  grid-generation  methods  In  use  today  can  be  grouped  Into  three 
categories:  those  which  apply  algebraic  shearing  or  stretching  of  the  coor¬ 
dinates,  those  which  are  generated  numerically  by  solving  a  Poisson  equation, 
and  those  which  are  baaed  on  conformal  mapping.  The  current  research  Is  con¬ 
cerned  with  an  example  from  the  third  category,  namely  the  transformation  Intro- 

3  9  10 

duced  by  Ives  and  Llutermoza.  ’  ’  Attempts  to  use  this  method,  as  reported  In 
Ref.  7  encountered  the  problem  of  non-orthogonallty,  especially  when  applied  to 
the  high  solidity  blade  rows  typical  of  present  designs  (slant  gap/chord  ratios 

on  the  order  of  0.7  or  less).  The  primary  goal  of  the  present  research  was  to 
extend  the  Ives-Liutermiza  technique  so  as  to  apply  to  turbine  blade  rows  of 
high  solidity. 


1 


As  the  first  step  in  achieving  this  goal,  an  Interim  Report^  was 
prepared,  documenting  a  number  of  improvements  that  had  been  made  since  issuance 
of  Ref.  9.  The  version  of  the  method  presented  in  Ref.  1  also  segregated  the 
portions  of  the  method  in  which  the  improvements,  to  achieve  orthogonal  grids 
at  high  solidity,  were  to  be  made. 

The  principal  modification  contained  in  the  present  report  is  the  re¬ 
placement  of  the  Theodorsen-Garrick  mapping  (used  as  one  of  the  steps  in  Refs. 
1,8-10)  with  a  derivative  form,  as  suggested  by  Bauer  et  al^^^^.  Use  of  this 
replacement  makes  it  possible  to  handle  blade  rows  of  high  solidity.  In 
addition,  a  further  step  involving  the  Schwarz-Christoffel  transformation  can 
be  used,  to  generate  a  fully  orthogonal  grid.  These  new  developments  are  de¬ 
scribed  in  Section  2  (along  with  an  updated  version  of  Ref.  1,  and  8-10,  for 
completeness) . 

The  potential  flow  through  a  cascade  of  blades  can  be  written  down 
algebraically,  once  the  cascade  has  been  mapped  into  a  unit  circle.  This  in¬ 
formation  has  been  added  as  Section  3,  using  one  of  the  mapping  steps  in  which 
the  blade-surface  image  is  mapped  into  a  circle. 

The  final  sections  contain  a  description  of  the  computer  code  and 
comments  on  its  range  of  applicability. 


Section  2 

CONFORMAL  TRANSFORMATION  METHOD 


The  method  of  Ives  §  Liutermoza  consists  of  a  sequence  of  transforma¬ 
tions,  which  map  a  two-dimensional  cascade  into  a  rectangle.  The  notation  and 
coordinate  system  used  to  define  the  cascade  are  shown  in  Fig.  1.  The  X- -co- 


The  quantities  5  ,  n.  and  the  angle  3^  denote  the  "streamwise,  normal”  coordinates, 
in  terms  of  which  the  blade  profiles  are  sometimes  defined.  These  reduce  to 
the  ^  set  if  W  is  taken  as  zero.  The  origins  of  both  of  these  coordinate 
systems  are  arbitrary. 

These  coordinates  define  complex  variables  3  and 

J?  »  S  -h  i  rt,  ^  =■  C2-1] 

The  points  ZN  and  ZT  are  taken  anywhere  near  the  centers  of  curvature  of  the 
leading  and  trailing  edges,  while  ZLE  and  ZTE  are  points  which  divide  the 
"pressure"  side  of  the  blade  (i.e.,  its  concave  surface)  from  the  "suction" 
side  (its  convex  surface) .  These  points  can  be  chosen  anywhere  on  the 


leading-  and  trailing-edge  contours;  ZTE  is  the  point  that  will  be  connected 
to  the  "point"  at  downstream  infinity  by  one  of  the  grid  lines,  if  that  option 
is  chosen  (i.e.,  ISHEAR=1). 

For  the  case  of  a  sharp  trailing  edge  (ITE=0),  ZT  must  equal  ZTE, 
and  for  a  sharp  leading  edge  (ILE=0),  ZN  must  equal  ZLE.  The  included  angle 
at  a  sharp  trailing  edge,  TT  ,  must  be  specified.  This  is  illustrated  in 
Fig.  2,  for  the  cascade  used  in  Ref.  7. 


Figure  2.  Blade  Geometry  of  Ref.  7 

This  blade  row  uses  the  NACA  65(12)10  profile,  has  a  slant-gap  chord  ratio 
S’C^/c  -  JL  ,  where  5Gj  is  the  slant  gap 
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and  a  stagger  angle  V  of  28.4®.  This  geometry  is  used,  below,  to  illustrate 
the  steps  in  the  Ives  transformation.  A  comparable  set  of  illustrations,  for 
a  cascade  of  turbine  blades  with  rounded  trailing  edges,  is  given  in  Ref.  10. 

Figure  2  shows  the  relation  between  the  trailing-edge  angle  Z  and  an 
exponent  EX  (which  is  used  in  one  of  the  transformation  steps  described  below). 
This  relation  applies  only  for  a  sharp  trailing  edge.  For  a  rounded  trailing 
edge>  is  not  related  to  the  180-degree  trailing-edge  angle,  but  is  chosen 
as  a  number  in  the  range  0.2  to  0.4,  as  described  below. 

The  blade  shape  is  input  as  two  tables  of  coordinate  pairs,  one  for 
the  pressure  surface,  and  one  for  the  suction  surface.  These  coordinates,  plus 
the  leading-  and  trailing-edge  points  ZLE  and  ZTE  are  then  arranged  in  an  array 
indexed  by  KJ,  where  1CJ=1  at  the  trailing  edge,  and  where  the  numbering  pro¬ 
ceeds  around  the  pressure  side  to  the  leading  edge  (KJLE),  and  then  along  the 
suction  side  to  the  trailing  edge,  where  the  point  denoted  by  KJMX  is  a  repeat 
of  that  denoted  by  KJ=1.  This  notation  is  shown  in  Fig.  3. 


ZLE  - 5 

(KJ  =  KJLE) 


ZTE 


- 


(KJ  =  ;  ^  KJ  MK) 


Figure  3.  Notation  for  Blade-Surface  Coordinates 

The  quantities  KJS  and  KJP  need  not  be  equal;  they  are  limited  to  a  maximum 
value  of  80  by  a  dimension  statement  in  the  current  version  of  the  program. 

The  first  step  in  the  Ives  transformation  is 
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The  fact  that  only  differences  of  i  -values  are  used  is  what  accounts  for  the 
arbitrariness  of  the  origin  in  Fig.  1.  On  the  suction  (  5  )  and  pressure  (  P  ) 
sides,  the  function  ^(2)  has  the  Fortran  equivalents: 
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The  arguments  of  the  sine  functions  can  be  written  in  a  simpler  form. 


as  follows: 
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where  ^2  stands  for  2-2^  or  Z -Z^  in  the  expressions  for  and  ,  respectively. 
By  noting  that  (see  Fig.  4) 
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it  follows  that 
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Similarly, 

‘  -‘  <>=-*«>}• 

Thus  the  function  ^  can  be  written  as 
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where 


w  ^t;  M 


The  shape  of  the  blade-surface  image  in  the  g-plane  can  best  be 
illustrated  by  considering  a  cascade  of  airfoils  consisting  of  flat  plates 
with  circular  leading  and  trailing  edges: 


(2-6) 


(2-7) 


FIGURE  5  NOTATION  FOR  FLAT- PLATE  CASCADE 


For  points  that  lie  on  the  circle  at  the  trailing  edge,  the  corresponding 

location  in  the  a  -plane,  to  first  order  in  ,  is 


The  variation  of  ^  with  V  and  is  shown  in  Fig.  6.  For  points 

thaft  lie  on  the  circle  at  the  leading  edge,  the  corresponding  relation,  to  first 

order  in  Q  ,  is 

1  TT 

Along  the  upper  surface  of  the  flat-plate  portion  of  the  airfoil,  the  coordinates 
in  the  physical  plane  have  the  values; 

C  ( '('  +  it’)  o  V 

^ 

where  T  is  measured  from  ^  from  ,  and 

r  'T  =  C 

The  resulting  formula  for  ^  shows  that  at  midchord  Ci-e-,  f  ~  ^  ^ 

These  general  features  of  the  ^-plane  are  shown  in  Fig.  7,  for  a 
negative  stagger  angle.  Note  that  a  clockwise  path  around  the  airfoil,  starting 
from  the  trailing  edge,  leads  to  a  clockwise  path  around  the  small  circle  in  the 
^  -plane,  followed  by  a  transit  to  the  large  circle  along  one  side  of  the  flat- 
plate  image,  and  then  a  counter-clockwise  path  around  the  large  circle,  and  back 


FIGURE  7  THE  MAPPING 


to  the  small  circle  along  the  other  side  of  the  flat-plate  image.  For  blades 
of  finite  thickness,  the  flat-plate  image  becomes  a  pair  of  lines,  which  fair 
smoothly  into  the  large  and  small  circles. 


In  writing  the  polar  angle  of  any  point  in  the  ^-plane,  care  must 
be  taken  to  retain  continuity  across  the  cut  (along  the  negative  real  axis) 
that  is  used  by  the  FORTRAN  ATAN2  function.  This  is  achieved  by  writing: 

I2£AL(<J)|  +  2  ~ 

where  the  'branch  number'  BR  is  set  equal  to  zero  at  the  trailing-edge  point 

(orkT^lfor  a  sharp  trailing  edge)  if  the  value  returned  by  ATAN2  at  this 
point  is  positive  (and is  set  equal  to  -1,  if  the  value  returned  is  negative) 
and  where  is  incremented/decremented  each  time  the  branch  cut  is  crossed 

in  the  counterclockwise/clockwise  direction.  For  a  sharp  leading  edge,  the 
increment  occurring  at  the  leading-edge  point  must  be  added  specifically,  since 
the  radius  of  the  large  circle  in  the  plane  is  infinite  for  this  case. 


The  next  step  is  the  exponentiation,  defined 
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Vk 
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(2-8) 


The  result  of  this  step  is  shown  in  Fig.  8,  for  the  cascade  depicted  in  Fig.  2. 
The  value  of  £%  for  a  sharp  trailing  edge  is  fixed  by  the  angle  ;  for  a 

round  trailing  edge,  a  value  of  in  the  neighborhood  of  0.2  to  0.5  will 

usually  reduce  the  variation  of  jfl  so  that  the  resulting  curve  in  the  _Q  ■ 
plane  will  have  a  polar  radius  that  is  a  single-valued  function  of  near 

the  trailing  edge. 


The  next  transformation  is  a  bilinear  one,  whose  purpose  is  to  produce 
a  curve  in  the  oJ  -plane  that  can  be  mapped  into  a  unit  circle  in  a  subsequent 


where  au ,  b  ,  and  C  are  complex  constants.  Three  conditions  must  be  assigned, 
in  order  to  evaluate  these  constants:  Ives  suggests,  for  two  of  them,  that 
the  images  of  upstream  and  downstream  infinity  be  mapped  into  cO"i  1  .  As 
the  third  condition,  he  recommends  that  the  centroid  of  the  blade-surface  image 
in  the  oO -plane  be  forced  to  be  close  to  the  origin.  This  condition  was  applied 
in  Ref.  9;  however,  it  has  been  found  simpler  to  impose  a  condition  on  the 
ratio  between  the  maximum  and  minimum  radii  in  the  aJ -plane,  as  outlined  below. 
The  calculation  of  the  centroidal  location  has  been  retained  in  the  present 
code,  for  informational  purposes.  (Details  on  how  this  is  calculated  can  be 
found  in  Ref.  9). 

The  locations  of  the  images  of  upstream  and  downstream  infinity  in 
the  ^  -  andil  -planes  can  be  expressed  as  follows:  in  general, 


As  >-*>-a>at  finite  (see  line(T)of  Fig.  9),  the  large-argument  approximations 
to  the  hyperbolic  functions  give 


Similarly,  as  y— ^♦■wat  finite  (see  line@of  Fig.  9),  the  corresponsing 
limit  is 


=  evlp  [  -  ^  (€r- 


These  are  written  as 


}C-K 


a  (- 


i 


In  terms  of  the  parameters  E  ,  F  ,  and  c  ,  the  transformation  can  be  written 


OJ  = 


-a  = 


f  -h  (Ec.-F'>  n, 

c  -  rL 

(  00  *  F  )  c  -  £ 
i*j  -  F  +  E  c. 


C2-15) 


(2-16) 


The  latter  relations  can  be  used  in  an  iteration  procedure  to  find  a  value  of 
C  that  will  minimize  the  ratio  RMAX/RMIN,  where  RMAX  and  RMIN  are  the  maximum 
and  minimum  values  of  10)1  over  the  set  defined  by  fCJ*l  to  KJMX.  The  iteration 
process  is  as  follows;  on  alternate  iterations,  values  of  c  are  chosen  that 
will  either  reduce  RMAX  or  increase  RMIN.  This  is  done  by  solving  Equation  2-16 
for  C  : 

XlcJ  -h  £  ~  F  XL 


oJ  4-  F  -  E  XL 


(2-17) 


On  the  first  iteration,  C  ^  1  is  used;  the  values  of  RMIN  and  the  index 
KJ=KJMN  at  which  it  occurs  are  then  found.  Next,  a  new  value  of  C  is  calculated, 
using  Eq.  2-17,  such  that  the  new  value  of  aJ (KJMN)  will  equal  1.1  times  the 
value  just  found.  For  this  new  mapping  into  the  oJ-plane,  the  value  of  RMAX 
and  the  index  KJ=KJMXX  at  which  it  occurs  are  found.  For  the  third  iteration, 
a  value  of  C  is  used  such  as  to  generate  a  new  value  of  cO  (KJMXX)  equal  to 
0.9  times  the  previous  one.  This  alternating  cycle  is  then  continued  until 
the  ratio  RMAX/RMIN  is  less  than  the  tolerance  RTOL.  This  tolerance  is  assigned 

a  default  value  of  3.0;  values  up  to  6.0  have  been  handled  successfully  by  the 

subsequent  steps  in  the  transformation. 

The  iteration  on  C  can  be  bypassed,  if  C  is  already  known,  by  setting 

IGOT=l  and  reading  in  the  value  of  c  •  Figure  10  is  the tO  -plane  for  the  cascade 

shown  in  Figure  2. 


•- 


-  J  .*  «**  .'->**  •’*  .* 
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Next,  the  blade-surface  image  in  the  -plane  is  mapped  into  the  unit 
circle  in  the  -plane  with  the  trailing  edge  at  ^  =  1  by  either  of  two  methods 
The  first  is  the  Theodorsen-Garrick  transformation: 


OJ 


L  ^ 


(2-1? 


To  determine  the  coefficients,  the  values  of  CJ  and  4  on  the  blade  surface 
are  written  as 


CO  =  r  C&)  e 


i  9 


,  ^  »  J  C 


(2-1! 


Then  the  real  and  imaginary  parts  of  the  transformation  are 


JLn  r 

^  r 

-  [ 

(p  -  3  4>~^ 

(2-21 

^  r  .  •  '  1 

6 

-  4>  + 

2;  \  PI-  j4>  -h  B:  ^  (p  \ 

(2-2 

FIGURE  12  INTRINSIC  COORDINATES 


The  coefficients  are  then  determined  by  the  following  iteration  procedure: 
an  equally-spaced  array  of  values  of  ^  is  set  up,  and  all  the  coefficients 
ar^'  initially  set  equal  to  zero.  Then  the  second  equation  gives  9=  cp  as 
the  first  approximation  for  9  .  For  each  of  these  values  of  ^  ,  a  corres¬ 
ponding  value  of  jir\  r  can  be  found,  from  a  spline  fit  to  the  coordinates  of 
the  blade-surface  image  in  theou-plane.  These  known  values  of  are  then 

used  in  the  first  of  the  equations  above,  to  find  the  next  approximation  to 
the  and  Sj  coefficients.  These  coefficients  can  then  be  used  to  give  the 
next  approximation  to  dC<p) ,  and  the  process  is  continued  until  convergence 

is  reached,  to  some  preassigned  tolerance.  Fast  Fourier  transform  tech- 
12 

techniques  can  be  used  in  the  processes  of  evaluating  the  second  equation,  for 
known  values  of  the  coefficients,  and  of  determining  the  coefficients  from 
the  first  equation  with  known  values  of  r  .  These  techniques  were  applied 

in  the  present  calculations.  The  details  are  given,  in  Appendix  A. 

Sufficient  conditions  for  convergence  of  this  iteration  process 

13 

have  been  discussed  by  Warschawski.  For  the  present  case,  these  conditions 
are  not  met;  in  particular,  it  is  required  that  the  maximum  and  minimum  values 
of  r  (9)  obey  the  relation: 

IRmAx  '  RMRX 

-  -  1  <  0.295  or  -  <  l.(o7e 

V  RmiN  RMlN 

This  condition  is  seldom  met  in  practice.  It  was  found,  however,  that  the 
iteration  process  would  converge  if.  a  relaxation  parameter  was  used.  i.e., 

12.  Cooley,  J.W.,  Lewis,  P.A.W.,  and  Welch,  P.D.,  "The  Fast  Fourier  Transform 
Algorithm:  Programming  Considerations  in  the  Calculations  of  Sine, 

Cosine  and  Laplace  Transforms",  Journal  of  Sound  and  Vibration  _12, 

(1970)  pp.  315-337. 

13.  Warschawski,  S.E.,  "On  Theodorsen's  Method  of  Conformal  Mapping  of 

Nearly  Circular  Regions",  Quarterly  of  Applied  Mathematics  3,  (1945) 


if  a  relaxation  parameter  was  used,  i.e.,  the  values  of  6  called  for  by  the 

second  equation  C  called  were  not  used  in  the  first  equation,  but  were 

replaced  by  B  ^  0. 1  6*^  *  O.S  Q  With  this  relaxation  factor,  the 

new  old 

iterations  were  convergent:  after  68  iterations,  the  maximum  change  in  any 
of  the  values  of  6  was  less  than  10  radians.  The  variation  of  &  with  <p  is 
shown  in  Fig.  11.  Calculations  for  other  cases,  not  shown  here,  have  required 
relaxation  factors  as  low  as  0.02  for  convergnece.  A  recent  review  paper  by 
Henrici  (Ref erence  14)  calls  attention  to  the  applicability  of  under-relaxation 
in  this  problem. 


When  the  solidity  is  high,  the  curve  of  9  vs  becomes  extremely  steep; 
typical  results,  for  a  gap/chord  ratio  of  0.8,  are  shown  in  Ref.  10.  For  such 


cases,  a  preferable  approach  is  to  use  a  derivative  form  of  the  transformation 


11 


To  determine  the  coefficients,  the  cJ  -plane  is  expressed  in  terms  of  the  intrinsic 
coordinates  V  and  /S>  ,  where  S  is  arclength  measured  from  the  trailing  edge, 
and  ^  is  the  angle  between  the  real  axis  and  the  tangent  to  the  curve  (see  Fig.  12) 
The  relations  between  the  intrinsic  coordinates  and  the  polar  angle  in  the 
-plane  are:  i 

-C  1-^1-  S  ! 


'(j 


re. 


since  (  and 


If  the  Fourier  coefficients  are  written  as 


;p  =  ^ 

j  * 


1 


Then  the  basic  equation  can  be  split  into  its  real  and  imaginary  parts  as: 


14. 


Henrici,  P.,  "Fast  Fourier  Methods  in  Computational  Complex  Analysis", 
SIAM  Review  21,  (1979)  pp.  481-527. 


= 


s  "DT 


<^  Ctfo  ^4>  +  T>l2i  A*— "1 4) 

r- 


The  coefficients  are  then  found  by  the  following  steps: 

1.  ^  and  S  are  found  from  the  numerical  data  that  define  the 
blade-surface  image  in  the  uJ  -plane 

2.  A  first  guess  at  /  is  made,  at  equally  spaced  points 

on  the  unit  circle  in  the  ^ -plane. 

3.  The  corresponding  values  of  the  arc length  are  found,  from 

^  -  f I  \  .lA 


r 


A  spline  fit  of  the  5  data  is  used  to  evaluate,  at  the 
equi-spaced  <f>  -values,  the  quantity: 


5.  The  Fast  Fourier  Transform  is  then  used  to  find  the  coefficients 
which  fit  the  data  of  step  4. 

6.  These  coefficients  are  then  used  [again  with  the  FFT)  to 
evaluate  the  next  approximation  to  jt^/ .  and  steps  3 
through  6  are  repeated  until  convergence  is  achieved. 

The  details  of  this  procedure  are  very  similar  to  those  involved  in 
the  Theodorsen-Garrick  technique;  the  FFT  steps  described  in  Appendix  A  can 
be  taken  over,  with  only  minor  changes. 

This  procedure  was  found  to  converge  in  15  iterations  (to  an  arclength 
error  less  than  .03)  without  relaxation.  The  resulting  variations  of  S'  and 
with  4  are  shown  in  Figure  13. 


The  next  transformation  is 


£-a 


(2-22 


where  0^-,  /3 ,  and  are  chosen  so  as  to  place  the  images  of  2  »  ±  <»  at 
=  ±  S  ,  while  the  blade  surface  continues  to  be  the  unit  circle .  These  images 
are  located,  respectively,  at  uJ»±t,  and  zj  -  Explicit  formulas 

for  cc  ,  /3 ,  and  in  terms  of  4 .  and  4.  are  given  by  Ives^  as 


VI  -*-  +  -J  I 


VIZ-  -vx  -  r 


^b)  -2/< 


~  ^  ^  *  Ab) 


*rtV-> 


Mokry  (Reference  151  has  pointed  out  that  the  formula  for  can  be  simplified, 
as  follows:  define 

C  5  - - 


s  '  ICI  -  Vlc/^-~r'  ,  y 

C  -  SICI  ^ 

OL  .  ^ - - - 3—  >  p 

c  -  SICJ  <8 


c  -  SICI 
idi 

Id  ~S c 

ic\  Aa  ~ 


(2-24) 


The  computer  program  described  below  does  not  use  these  simplifications. 

Next,  it  is  necessary  to  find  and  ,  given  uJ  ■  -  1.  This 
was  done  by  Newton- Raphson  iteration: 


_G££2. 

dG 

cLi, 


where  for  the  Theodorsen-Garrick  procedure 


SCO  » 

i, 

ct  Gr 

w  r 

oL 

^  [ 

15.  Mokry,  M.,  "Comment  on  Analysis  of  Transonic  Cascade  Flow  Using  Con¬ 
formal  Mapping  and  Relaxation  Techniques",  AIAA  Journal  16,  No.  1, 
(January  1978)  p.  96. 


In  order  to  find  an  initial  guess  for  2;  ,  a  preliminary  calculation  was  made, 
for 

1 2;  I  »  0.49  (.1)  0.99,  arg?  »  -60“  flO")  ♦  60“ 

From  this  set,  the  value  of  25  which  gave  to  nearest  to  +1  was  chosen  as  the 
initial  guess.  This  set  was  then  repeated  with  ?  replaced  by  -  ^  ,  to  obtain 
the  value  of  5  that  gave  a)  nearest  to  -1.  The  locations  of  key  points  in  the 

4  and  planes  are  shown  in  Fig.  14. 

When  using  the  derivative  form  of  the  cJ^  ^  transformation,  it  is 
necessary  to  integrate  the  derivative  ;  this  was  done  using  Simpson’s 

Rule,  starting  at  the  trailing  edge: 

I 

The  equations  used  in  this  case  for  the  Newton- Raphson  procedure  are 


V 

i 

I 


When  the  values  of  2^^  and  are  known,  the  blade-surface- image  points 
can  be  mapped  from  the  ?  -plane  to  the  -plane.  In  both  planes,  these  points 
are  located  on  unit  circles,  so  it  is  only  necessary  to  interpalate  in  Fig.  12 
to  find  the  ^  -values  for  the  points  defining  the  blade  surface.  This  is  done 
with  a  call  to  the  spline-fit  subroutine,  and  the  values  returned  by  it  are 
replaced  by  linear  interpolation  in  regions  where  the  (9- ,  <P  curve  is  so  steep 
that  the  spline  fit  returns  non-monotonic  values. 

The  final  transformation  now  uses  an  elliptic  function  to  map  the 

unit  circle  in  the  -plane  (with  cuts  along  the  real  axis  from  ±5  to  the 
circle)  into  a  rectangle: 

^  S  Sn  C  ^  ^  -k  ) 

where  the  parameter  ^  is  given  by 


(2-25) 


4  "  PlR/\)£ 


q  -  PLF)KJ£ 


Figure  14  The  ^  and  ^  Planes 


The  inverse  of  this  transformation  is 


V5 


cit 


V?  -  t*‘  Vj  -  ^  (2-27) 

This  latter  transformation  is  used  to  find  the  images,  in  the  4  -plane,  of 
the  blade-surface  points.  This  requires  an  expression  for  the  real  and  imag¬ 


inary  parts  of  the  incomplete  elliptic  integral  of  the  first  kind;  convenient 
formulae  for  this  purpose  are  given  by  Nielsen  and  Perkins^^  who  show  that. 


if 


S  CZ  ^  I  S) 


(2-28) 


then 

i  ) 

(2-29) 

where 

V;  ~S*’ 

(2-30) 

16.  Nielsen,  J.N.  and  Perkins,  E.W.,  "Charts  for  the  Conical  Part  of  the 
Downwash  Field  of  Swept  Wings  at  Supersonic  Speeds",  NACA  Technical 
Note  1780,  (December  1948),  Appendix  C. 


cr  «  [r*  + <S*- +  X  +  j  -  aa* 


(2-31) 


These  formulae  are  equivalent  to  those  following  Eqn.  115.01  of  Byrd  and 
Friedman;  the  formula  for  X  has  heen  rearranged  slightly  from  the  form 
given  by  Nielsen  and  Perkins,  to  avoid  the  occurence  of  negative  values 
under  the  square  root  sign,  which  can  sometimes  happen  in  the  numerical  evalu¬ 
ations  when  6  0  .  These  formulas  are  correct  along  the  branch  cuts, 

where  6*0  and  |q|/<S'  >  1  . 


and  where  F  denotes  the  incomplete  elliptic  integral  of  the  first 

kind,  with  real  arguments  A  and  given  as  functions  of  T ,  d  and  4,  : 


id 


X  -  [/  tr’*  S‘-  Vo  -r’/  *  cf‘(f  *♦  2  ♦  2T*/  ]  + 

-  VCr  -  ■,  *W]  /  H‘r‘ 


-  W  w  . 

S'VtV-', 


Numerical  evaluations  of  these  elliptic  integrals  were  done  using 

18 

twelve  terms  in  the  formulae  of  Luke: 


« 

s  J 


-t^  Vi  -yfc*  t' 


1  -  4^ 


—  \4> 

25  L 


+  22: 


(T^  4>y 


(2-32) 


17.  Byrd,  P.F.,  and  Friedman,  M.D.,  Handbook  of  Elliptic  Integrals  for 
Engineers  and  Physicists.  Springer  Ver lag,  Berlin  (1954). 

18.  Luke,  Y.L.,  "Approximations  for  Elliptic  Integrals",  Mathematics  of 

Computation,  ^  (1968)  627-634. 


CvlUll 
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where 

<t>  <  rt/2  ,  V;  ,  (9^  ^  ^  rr/55  • 


I 


For  the  case  where  <p  =  ri/2  (the  complete  integral)  the  approximate  formula  is 

‘  -r]  (2-33: 

1-  KYI  S  ♦  YY\ 

The  signs  in  the  formulas  above  pertain  to  the  first  quadrant  in 
the  f\  -plane  i  Z  O  ,  6  i  O ) ,  which  maps  into  a  rectangle  in  the  first  quad- 
rant  of  the  ^  -plane,  with  sides  located  on  the  lines 


K(4) 


oLt 

- 1* '  vrrp^T^ 


Y  kVi) 


_ _ 

V;  -  t*‘ 


(2-34) 


The  remaining  three  quadrants  in  the  /]  -plane  map  into  the  remain- 
ing  three  quadrants  in  the  ^  -plane,  as  described  in  Reference  (19),  p.  377;  the 
cuts  from  t  S  to  the  unit  circle  along  the  real  axis  become  the  left  and 
right  sides  of  the  rectangle  in  the  ^  -plane: 


Erdelyi,  A,,  et  al.  Higher  Transcendental  Functions,  Volume  2,  p.  377, 
McGraw-Hill  Book  Company,  New  York  (1953) . 
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Kj  -  PuflNB 


^  -  PLRNE 


Figure  15  The  and  ^  Planes 


Finally,  the  ^  -  plane  is  re-normaliied,  so  as  to  lie  between  -1  and  +1  on 
both  axes: 


K(k)  ^  ^  K'L-k) 


(2-35 


A  grid  is  now  to  be  set  up  in  the  ^  plane,  and  mapped  back  to 
the  Z  -  plane.  This  process  is  facilitated  by  first  rearranging  the 
quadrants  in  the  ^  -  plane,  by  using  the  periodicity  of  the  elliptic  functions 
as  follows:  in  the  first  and  second  quadrants,  let  ^  ^  ,  and  in  the  third 

At  ^ 

and  fourth  let  ^  ~  -  f  Z  K  (A)  )  and  use  the  relations  (see  for 
example,  Eqs.  122.00  and  122.04  of  Reference  17. 


-'T, 


5,,  «-SWCf  =  Srt[-Cf+2K)]  =  Sn^-K-(^+K)j  (2-36) 


The  C  plane  then  has  the  form: 


Kc^) 


■3Ka) 


-2X'(*) 


K-tJkJ 


- 1  o. 


Figure  16  The  ^  Plane 


Two  types  of  grid  can  be  selected  in  the  ^  plane:  a  rectangular  one 
(if  ISHEARaO)  or  a  sheared  one  (if  ISHEAR=1).  The  latter  is  the  default. 


For  the  rectangular  grid,  equally  spaced  points  are  assigned,  ac¬ 


cording  to 


,  K-  (2.37) 

~  i  XU)-(L-0  L  .  ,I.MX 


where 


^K(-k) 
KMX  -? 


>  ^  K  it 


k  K\-k) 

i  '  L/^X-1 


(2-38) 


Becuase  the  mapping  is  conformal,  the  images  of  these  grid  lines  will  intersect 
at  right  angles  when  mapped  back  to  the  physical  plane. 

The  rectangular  grid  has  the  property  that  in  general  the  trailing 
edge  is  not  connected,  by  a  grid  line,  to  the  point  at  downstream  infinity. 
However,  such  a  connection  is  a  desirable  feature  in  certain  flowfield  codes 
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'  -  •  -  -  - • .A  -  -  - 


O  *-  • 


»  i 
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'  *►  v 
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(see  Ref.  7,  for  example).  To  allow  this  feature,  the  ISHEAR=1  option 
establishes  a  sheared  grid,  consisting  of  the  same  lines  as  above,  but 

A 

replacing  the  lines  by  a  set  of  parabolas  which  intersect  the  blade 

surface  at  90  degrees,  and  are  displaced  from  a  base  parabola  that  connects 
the  trailing  edge  to  the  image  of  downstream  infinity: 


This  grid  is  given  by 


iK'<*  > 

?,(TE)  .  (X-,) 

T , 


,  A  - 


L^X 


Const 


K  •  ;,2, 
L  »  t,2, 


,  X/^X 
^  L  MX 


where 


const 


-  [k  '(-k)] 

^  [3  K  (-Jt)  +  <T£)2 


^/TE)  -  fU  {^[kJ  -  l]Y 


Some  points  on  these  parabolas  will  lie  outside  the  range 

-  3K  (*)  ^  ^  +  K(-^) 

When  this  occurs,  equivalent  points  are  found  by  adding  or  subtracting 
the  period  of  the  elliptic  sine.  In  addition,  the  base  parabola  is  always 
joined  to  the  lower-left  corner  of  Fig.  15,  by  subtracting  4/<(^Jfrom  the 
real  part  of  ^  ,  if  the  latter  is  greater  than 

An  alternate  method  01  joing  the  images  of  the  trailing  edge  and  down¬ 
stream  infinity,  while  retaining  an  orthogonal  grid,  is  described  in  Section  4. 
It  involves  use  of  the  Schwart-Christoffel  transformation. 


This  completes  the  definition  of  the  grid  in  the  ^  -plane.  Each 
of  these  grid  points  must  now  be  mapped  back  to  the  physical  plane.  The 


first  transformation  is: 


q  =  ^sn(^  i  ;  «-) 


C2-41) 


The  elliptic  sine  of  a  complex  argument  is  expressed  as  (Ref.  11,  Eq.  125.01) 

5n  (UL  +  LV, 

5n  Cu-.-t)  dLndr,  -li')  +  i  CnCu-,-k)  cLnCa.,*)  sn  (  v.-ii')  c.nCV,-k') 

T  -  Sr)^  Cv,-k")  aLn^iU-.-k)  C2-42) 

The  functions  in  this  expression  are  evaluated  by  the  Arithmetic-Geometric 
Mean  method  (see  Ref.  20,  p.  571): 

Set 


i  i 
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..  •  •  •  1.- 
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«.  *  V  •  • 

%v-  •. 

^“,1  -*4 

>;-v 
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and  then  calculate 


2  ^^n-t  ^n-f)  > 


'a.  b 

n  .  f  n-  1 


~  -  C^n-i~  1 ) 


(2-43) 


until  Cn  =  0  to  a  prescribed  tolerance  (io~^  '^^.s  used  in  the  present  case.) 


Then  form 


Q  ^  -  2  •  <X„  •  u- 


and  calculate  >■■  ) 


Pn-1  =  J  ^ 


(2-44) 


Then  the  desired  results  are  given  by 


« 


20.  Abramowitz,  M.  and  Stegun,  I. A.,  Handbook  of  Mathematical  Functions, 
National  Bureau  of  Standards,  Applied  Mathematics  Series  55  (1964). 
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3n(u.-lt)  ~  ^  (p^  ^  cnCct.Tfe)  =  c^  <p^ 


oLn  (u.,-Jk)  =  _ ^ 

CM  id},  -  4). ; 


(2-45) 


.‘•V.-.N-C* 

v'-I-'-j-' 

*.  *1  _  w*„  «•, 


These  values  of  r]  are  then  mapped  back  to  the  ^  -  plane 
by 

^  n-  3' 

and  thence  to  the  cu-  plane  by  either 


OJ  =  ^  ^  2  ( ^  j  *  i  Sj  )  ^  1 

^  d  »o  *  J 


(2-46) 


(2-47) 
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if  the  Theodorsen-Garrick  transformation  is  benig  used  or  by 

.r  y. 


oJ  ^  ^ 


‘sr/»«.T' 


if  the  derivative  form  of  the  cJ^  \  transformation  is  being  used.  The  values 
used  for  the  starting  point  of  this  integration  are  as  follows:  at  first,  the 
equation  is  integrated  frora<o--(  to  -  + 1  ,  (,  ^  =  ’JgVo^^^so  as  to  establish 
the  periodic  boundary.  Then  a  series  of  integrations  is  done,  at  constant  K  , 
for  L  ranging  from  the  blade  surface  to  the  periodic  boundary  (the  values  of 
?  and  iJ  on  the  blade-surface  image  were  found  earlier  in  the  FFT  procedure). 

These  integrations  are  overdetermined,  in  the  sense  that  both  end  points 
are  known,  as  well  as  the- values  of  the  derivatives  in  between.  Moreover,  the 
integrations  used  to  establish  the  end  points  have  been  done  by  several  techniques 
(among  them,  Simpson's  rule  and  the  FFT  procedure).  Thus  the  integrations  from 
a  given  starting  point  do  not  always  terminate  at  the  desired  end  point;  when 
this  occurs,  a  small  adjustment  is  made,  in  the  values  of  the  derivatives,  so 
as  to  guarantee  the  proper  end  points. 


.'rV.V'.'ji 


V 
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Finally,  the  value  of  is  found  by  inverting; 


whose  inverse  is 


The  final  result  of  this  process  is  shown  in  Figure  17,  for  the  cascade  of 
Figure  2.  This  grid  was  found  using  the  derivative  form  of  the  transformation; 
a  comparable  grid  found  with  the  Theodorsen-Garrick  technique  is  shown  in 
Ref.  7. 

For  higher  solidities,  the  derivative  form  of  the  "S  transformation 
can  be  used;  Figure  18  shows  the  grid  that  results  for  a  gap/chord  ratio  of 
0.5.  The  grid  points  calculated  for  the  region  near  the  leading  edge  have  been 
omitted  from  this  figure,  because  they  cannot  be  located  with  sufficient  accuracy. 
The  basic  reason  for  this  loss  of  accuracy  can  already  be  seen  at  a  gap/chord 
ratio  as  high  as  1.0  (see  Fig.  13):  the  portion  of  the  blade-surface  image 
near  the  leading  edge  has  a  nearly  vertical  slope  in  this  figure.  Thus  while 
the  derivative  form  of  the  transformation  does  coverage  for  high  solidities, 
the  grid  which  it  leads  to  suffers  from  a  loss  of  accuracy  near  the  leading 
edge. 

This  loss  is  aggravated,  in  the  present  case,  by  the  adjustment  of 
derivatives  that  is  made,  to  enforce  periodicity.  Near  the  leading  edge,  the 
adjustments  required  become  significant,  and  their  cumulative  effect  is  to  pro¬ 
duce  a  non-orthogonality  at  the  blade  surface.  Further  study  is  required,  in 
order  to  minimize  this  problem  of  inaccuracy  near  the  leading  edge. 


wSSBA 


POTENTIAL  FLOW  FIELD 


It  is  possible  to  write  down  in  algebraic  form  the  complex  potential 
for  inviscid,  incompressible  flow  in  the  -plane,  and  this  solution  then  gives 
the  flow  in  the  physical  plane.  The  complex  potential  is  (see,  for  example. 


Ref.  21). 

W  ■=■  4>  f  (1  'f' 
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where  (p  and  «.j>  are  the  velocity  potential  and  stream  function,  CJ  and  are 
source  and  vortex  strengths,  oC.  is  the  angle  of  attack,  and  S  the  parameter 
described  earlier.  The  complex  velocity  is; 


The  Kutta  condition,  requiring  that  at  gives  the  following 

relation  for  ^  ^  Q  ; 


Putting  this  back  into  the  expression  for  /  JL'A.  ^  is  possible  to  find 
a  second  stagnation  point  (near  the  leading  edge)  at 


The  values  of  ^  and  Q  are  now  found  from  conditions  on  the  velocity 
at  infinity;  in  general,  the  velocity  components  U.  and  u'  in  the  plane 

are  given  by: 
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,  +<^>  the  image  in  the  i^^-place  is  T^^-t-  S;  then: 
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By  using  the  approximate  form  of  the  Jl^  i^y  relation  as  -*  a?  ^  it  is 
possible  to  show  that 


_n.  -  JL  (.00) 
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A  similar  analysis  shows  that 

(“-“'l-a,  =  ■'I-] 

To  satisfy  the  continuity  requirement  that  U^-co  —  ^-ao  take  Q.  = 


cr  sc,/ 

fit 


Then: 


U-ao  -  ‘^  +  ce  =  d:  CiTQ.  of. 


+.  CC  [— '^e&“  ^l.~  ^  ~  ^/<5  ^ 


These  relations  are  shovm  in  Figure  19;  note  that  o4  is  <■  O  for  positive 
incidence. 

Once  the  constants  and  Q  are  established,  the  velocities  at  any 
point  in  the  flow  can  be  found,  by  using 

^  ^  M-Ti. 

d-Z  d-Ci.  J. 

Each  of  the  derivatives  in  this  formula  can  be  evaluated  exactly.  A  typical 
result  is  the  streamline  pattern  shown  in  Figure  20. 

This  solution  is  of  value  in  its  own  right  (especially  since  it  is 
found  from  purely  algebraic  formulae),  and  can  also  be  useful  as  a  set  of 
initial  conditions  for  a  time-marching  or  iterative  solution  of  the  flow 
pattern. 


Section  4 
ORTHOGONAL  GRIDS 


It  was  pointed  out  above  that  many  current  flowfield  codes  require  a 
grid,  in  the  computational  plane,  having  the  property  that  the  images  of  the 
trailing  edge  and  the  point  at  downstream  infinity  lie  at  comers  of  the  grid. 

This  property  is  not  realized  by  the  Ives  transformation;  the  point  at  down- 
stream  infinity  occurs  at  a  comer  of  the  ^  -plane,  but  the  trailing  edge 
does  not.  The  sheared  grid  discussed  in  Section  2  is  one  way  to  achieve  the 
desired  property,  but  it  has  the  disadvantage  that  the  resulting  grid  is  not 
orthogonal . 

The  Schwarz-Christoffel  transformation  can  be  used,  in  place  of  the 

(22-25) 

shearing,  to  achieve  an  orthogonal  grid.  A  number  of  recent  publications  ^ 

have  discussed  the  use  of  this  transformation;  its  great  versatility  offers 
several  possible  ways  to  use  it,  within  the  Ives  Transformation.  For  example, 
it  would  be  possible  to  replace  the  elliptic-integral  step  between  the  and 
^  planes  by  a  hyper-elliptic  integral  (which  is  a  form  of  the  Schwarz-Christoffel 
formula).  However,  a  much  more  straightforward  application  has  been  made  here, 
namely  to  map  the  parallelogram  formed  by  joining  the  trailing-edge  and  down- 
Stream  infinity  points  in  the  ^  -plane  into  a  rectangle.  This  step  can  be 
used  on  other  grid  generators  having  similar  properties,  such  as  that  of  Ref.  26, 
for  example. 

Figure  21  shows  the  parallelogram  formed  by  joining  the  images  of  the 
trailing  edge  and  the  point  at  downstream  infinity  in  the  ^  plane: 


The  particular  values  of  the  coordinates  shown  on  this  figure  are  those  peculiar 
to  the  Ives  transformation;  the  method  applies  to  arbitrary  locations. 


This  figure  is  mapped  into  the  upper  half  of  the  "t  -plane  by  the 
Schwarz-Christoff el  transformation : 
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where  the  parameters 


through  are  pure  real : 


The  exponents  in  the  above  equation  guarantee  that  an  integration  path  along 
the  real  axis  in  the  i  -plane  will  produce  a  four-sided  figure  with  the  proper 
angles.  The  four  paraimeters  through  and  the  complex  constant  M 

must  be  chosen  so  as  to  produce  the  desined  parallelogram. 


It  turns  out  that  only  two  parameters  are  needed:  can  be  set 
equal  to— tp  ,  to  —  ‘t^  ,  and  can  be  set  equal  to  1.0 
without  loss  of  generality: 
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This  set  of  parameters  produces  a  genuine  parallelogram,  i.e.,  the  lengths  of 
opposite  sides  are  equal.  Thus,  all  that  is  needed  is  to  select  so  as 
to  give  the  desired  ratio  of  the  adjacent  sides,  and  then  set  the  parameter  M 
so  as  to  give  the  proper  absolute  size. 


Numerical  integrations  of  this  formula  can  be  carried  out  using  the 
method  of  Ref.  22;  the  integral  over  a  distance  can  be  written  as: 
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and  where  and  P  are  the  parameters ; 
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The  length/width  ratios  of  the  parallelograms  produced  by  typical  cascades  are 
generally  5.0  or  larger.  Numerical  work  using  the  formula  above  shows  that 
such  ratios  require  values  of  that  are  only  slightly  larger  than  1: 

Asymptotic  formulas  for  the  lengths  of  the  slant  side  (called  side  1)  and  the 
horizontal  side  (called  side  2)  of  the  parallelogram  can  be  found  by  taking 
the  limit  as  €—^o  : 

"  ‘  -Av-  - '  +  AAt 

—  A?')  =  f  (t+l<-6)  (i*') 

The  change  of  variable  leads  to: 

where 
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The  leading  term  can  be  expressed  in  terms  of  gamma  functions;  using  the 
relations  (Ref  27^ 
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leads  to 
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To  find  the  analogous  result  for  side  2: 
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In  the  first  of  these  integrals,  let  t  +  ,  and  in  the  second,  let 

;  this  leads  to 

Vf  -i^-Z/tt  -/3/^  r  /  „ 

^  J  ^  [  1 

O 

The  first  of  these  integrals  can  be  estimated  by  an  integration  by  parts, 
while  the  second  can  be  approximated,  for  6-  -a-o  ,  as  a  Mellin  transform; 
the  leading  term  in  each  integral  is  jUK.6  ,  with  the  result: 


I 
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The  practical  application  of  these  asymptotic  formulas  is  to  provide 
a  first  guess  at  €  for  a  given  ratio  of  the  length  of  the  adjacent  sides  of 
the  parallelogram.  Values  of  6  on  the  order  of  10  ^  are  needed,  for  side 
ratios  of  5  or  so;  these  small  values  required  comparably  small  step  sizes 
in  the  numerical  evaluations,  in  order  to  preserve  accuracy. 


Once  the  parameters  for  mapping  of  the  parallelogram  are  established, 
the  next  step  is  make  the  upper  half  of  the  i:  -plane  into  a  rectangle;  this 
is  achieved  by  using  the  same  values  of  ^  and  M  with  =  1/2. 

An  orthogonal  grid  can  then  be  set  up  in  the  rectangle  plane,  and 
mapped  (numerically)  to  the  -t  -plane.  This  grid  can  then  be  mapped,  again 
numerically,  from  the  i  -plane to  the  parallelogram  in  the  ^  -plane,  and 
then  back  to  the  physical  plane,  as  explained  in  Section  2  . 


Section  5 
METRIC  EVALUATIONS 

The  coordinate  transformation  enters  the  flowfield  solution  algorithm 
only  in  the  metric  derivatives.  These  can  be  evaluated  by  differencing  the 
coordinates  themselves,  or  in  the  case  of  a  conformal  transformation,  by  evalu¬ 
ating  the  analytic  expressions  for  them.  These  analytic  expressions  are  derived 
in  Ref.  1,  and  the  code  listed  in  that  report  contains  the  Fortran  statements 
required  to  evaluate  these  expressions.  However,  as  pointed  out  in  Ref.  7^ 
the  truncation  error  resulting  from  the  use  of  analytic  metrics  in  the  finite- 
difference  flowfield  code  is  large  enough  to  cause  major  instabilities  in  the 
solution  algorithm.  It  is  highly  preferable  to  use  metrics  that  are  found  by 
differencing  the  coordinates  in  the  same  manner  as  the  flowfield  variables  are 
differenced.  A  progreim  to  achieve  this,  for  the  grid  conventions  used  in 
Ref.  7,  is  given  in  Appendix  D. 


Section  6 
COMPUTER  PROGRAM 


A  listing  of  the  computer  program  is  given  in  Appendix  B,  and  a 
dictionary  of  variables  is  given  in  Appendix  C.  This  section  contains  a 
general  description  of  the  program,  plus  some  specific  details. 

In  order  to  handle  complex  arithmetric,  all  variables  beginning 
with  the  letter  Z  are  declared  to  be  complex  by  an  implicit  type  specification 
at  the  beginning  of  the  program. 

The  input  is  generally  described  by  comment  cards  in  the  deck.  Cer¬ 
tain  blade-shape  parameters  must  be  read  in:  EX,  G,  H,  ZLE,  ZTE,  ZN,  ZT.  Also, 
if  IG0T=1,  ZC  must  be  read  in.  The  blade  shape  itself  is  defined  by  pairs  of 
coordinates  in  the  S  ,  n  plane  -  fCJS  on  the  suction  side  and  KJP  on  the  pressure 
side.  In  the  version  shown  here,  these  were  read  in  as  a  table  in  subroutine 
SHAPE. 

The  blade  surface  coordinates  are  numbered  from  1  to  KJMX;  point  number 
1  is  the  trailing  edge,  2  through  KJIM  are  on  the  pressure  side  from  trailing 
edge  to  leading  edge,  point  KJLE  is  the  leading-edge  point  (ZLE),  KJLP  through 
KJMXM  are  on  the  suction  side  from  leading  edge  to  trailing  edge,  and  point 
KJMX  repeats  the  trailing  edge. 

The  complex  sine  functions  are  calculated  next;  then  the  iterations 
to  determine  the  parameter  C  are  done.  The  initial  guess  provided  for  C  is 
Zc  =  1  -f- 1  ,  It  may  happen  in  some  cases  that  a  better  guess  is  required: 
in  particular,  it  is  necessary  that  the  value  of C  must  lie  outside  the  blade- 
surface  curve  in  the.Il  -plane.  If  it  does  not,  then  the  interior  of  this 
curve  in  the  A  -plane  is  mapped  to  the  exterior  of  the  blade-surface  image  in 
the  tO -plane.  This  fact  can  be  seen  from  the  discussion  by  Kober  (Ref.  28 
of  the  bilinear  transformation  applied  to  circles. 


28.  Kober,  H.,  Dictionary  of  Conformal  Transformations,  Dover  Publications, 
New  York  (1957) . 


After  the  parameter  C  has  been  found,  the  transformation  to  the 
plane  is  carried  out,  using  the  fast  Fourier  transfrorm  procedure  (see 
Appendix  A] .  The  actual  calculations  of  the  Fourier  coefficients  are  done 
in  subroutine  FFT2,  a  proprietary  program  of  International  Mathematical  and 

Statistical  Libraries,  Inc.  (IMSL).  This  routine  computes  the  fast  Fourier 
transform  of  a  complex  vector  of  length  equal  to  a  power  of  two  (here  2^) . 

The  coefficients  of  the  input  vector  are  given  in  normal  order  by  the 
array  named  as  the  first  argument  of  the  call;  the  coefficients  of  the 
output  vector  are  overstored  in  this  array,  in  reverse  binary  order.- 
The  subroutine  SHUFL  is  then  used  to  restore  this  output  to  the  normal  order. 
The  coefficients  in  the  series  expression  for  ^  are  determined  iteratively 
in  a  relaxation  process  that  is  terminated  when  the  maximum  change  in  O'  falls 
below  the  tolerance  ANGERR,  or  when  IMX  iterations  are  done.  If  the  derivative 
form  of  the  transformation  is  being  used,  the  maximum  change  in  S  is  required 
to  be  less  than  ANGERR. 


In  doing  these  iterations,  it  is  necessary  to  know  values  of 
at  given  values  of  &  ;  these  are  found  by  a  spline  fit  in  subroutine 

CISPLN,  which  is  a  straightforward  implementation  of  the  formulas  given  by 
29 

Ahlberg,  et  al.  Similarly,^  at  eiven  value.<!  of  s  is  found  by  a  spline 
fit. 


In  certain  cases  (typically  when  RMAX/RMIN  is  large)  the  calculated 
variation  of  d  with  may  be  non-monotonic;  if  this  occurs,  the  calculation 
should  be  repeated,  with  a  smaller  relaxation  factor.  The  progress  of  the 
iterations  is  printed,  showing  at  each  iteration  the  largest  change  in 
0  and  the  number  of  reversals  (i.e.,  the  number  of  occurrences  of  non-monotonic 
variation).  The  iterations  using  the  derivative  form  usually  do  not  display  a 
similar  problem,  and  converge  very  quickly,  with  OM  =  1. 


29.  Ahlberg,  J.H.,  Nilson,  E.N.,  and  Walsh,  J.L.,  The  Theppr  of  Splines 
and  Their  Applications  Academic  Press,  New  York  (1967) . 
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Next,  the  parameters  and  are  found,  starting  with  "best-guess" 
values  calculated  in  the  sectors  described  in  Section  2.  Once  these  are 


found,  the  mapping  to  the  -plane  follows.  The  calculations  that  link  the 
-plane  and  the  ^  -plane  are  done  in  subroutine  OMETA,  which  sums  the 
Theodorsen-Garrick  series,  using  complex  arithmetric.  This  subroutine  has 
been  modified  slightly  from  that  appearing  in  Ref.  9,  as  follows:  the  previous 
code  evaluated  sums  of  the  form 


US 


SUM 


z  2CC(JP)  t; 
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(4-1) 


by  a  sequence  of  multiplications  and  additions: 


ISUM  =  2CC  (US)  ^  ZccCi.^) 

The  current  version  uses  (Ref.  30,  p.  28) 


+  2cc(T)  (4-2) 
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Finally,  the  blade-surface  image  is  mapped  into  the  ^  plane,  using 
the  elliptic-function  fomnulas  of  Nielsen  and  Perkins^^  and  of  Luke,^®  as 
outlined  in  Section  2.  This  completes  the  mapping  of  the  blade  surface. 


It  is  now  possible  to  set  up  a  grid  in  the  ^  -  plane,  and  map  it 
back  into  the  physical  plane.  This  involves  straightforward  evaluations  of 
the  transformation  functions.  The  only  complication  is  the  need  to 
evaulate  the  Jacobian  elliptic  sine  (done  in  subroutine  JCELFN), 


30  Hartree,  D.R.,  Numerical  Analysis,  2nd  Edition,  Oxford,  Clarendon 
Press  (1958). 
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The  calculation  of  the  images  of  the  grid  points  in  the  various 
planes  is  bypassed  for  the  points  at  upstream  and  downstream  infinity,  and  at 
the  trailing  edge.  (The  trailing-edge  point  will  be  a  grid  point  if  I£HEAR=1). 
The  4  -plane  locations  of  upstream  and  downstream  infinity  are  arbitrarily 
assigned  to  finite  locations  given  by  linear  extrapolation  from  the  two  adjacent 
L  -values. 

The  point  at  upstream  infinity  will  be  a  grid  point  only  if  KMX  is 
odd;  in  this  case  I0E=1,  and  the  image  calculations  are  bypassed. 

Adjustments  to  the  grid-point  image  locations  in  the  Z  -plane  are 
sometimes  required  for  points  on  the  periodic  ouundary  (L=LMX)  for  values  of  K 
near  1,  KMX/ 2,  and  KMX.  At  these  points,  the  formula  used  in  going  from  the 
-CL  -plane  to  the  4  -plane  sometimes  cannot  distinguish  between  points  that 
are  separated  by  H  + t  Q  .  The  problem 'can  be  seen  best  in  the  cJ  -plane. 

(The  sketch  below  is  for  an  even  value  of  KMX;  the  same  picture  applies  for  an 
odd  value) : 


Note  that  the  same  point  in  the  oJ  -plane  (and  thus  also  in  the  il -plane)  can 
map  into  either  of  two  points  in  the  £  -plane,  which  differ  by  H  +  i&  .  The 
selection  is  guided  toward  the  correct  value  by  starting  on  the  blade  (L=l) 
and  working  toward  the  periodic  boundary  (L=LMX),  but  it  can  happen  that  the 
wrong  branch  is  chosen  during  the  iterations.  To  avoid  this,  the  imaginary 
parts  of  2  and  EA/  are  compared,  for  K  values  near  the  leading  edge,  and  the 
imaginary  parts  of  H  and  ZT  are  compared  near  the  trailing  edge,  and  the 
quantity  L‘S&  is  added  or  subtracted  (depending  on  the  value  of  K  )  where 
necessary. 

Finally,  the  real  and  imaginary  parts  of  £  are  written  to  unit  7 
(if  PNCHZA=TRUE) .  These  values  can  be  used,  in  a  separate  program,  to  calculate 
the  metrics  of  the  transformation  (see  Appendix  D) . 


53 


Section  7 

CONCLUDING  REMARKS 


The  method  described  above  is  capable  of  generating  grids,  and  the 
associated  incompressible,  inviscid  flow-field,  for  blade  rows  with  gap/chord 
ratios  as  low  as  0.8.  Below  that  value,  the  Theodorsen-Garrick  mapping  may  not 
converge.  In  such  cases,  a  derivative  form  of  the  transformation  usually  will 
converge,  but  the  grid  that  results  may  not  be  useful,  in  the  region  near  the 
leading  edge.  Further  development  of  this  part  of  the  technique  is  required. 

As  a  separate  mapping  step,  the  Schwarz-Christoffel  technique  has  been 
applied  to  the  problem  of  generating  orthogonal  grids.  Based  on  the  analysis 
given  above,  several  grid-generated  techniques  can  now  be  generalized  to  give 
orthogonal  grids. 
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APPENDIX  A 

DETAILS  OF  THE  FAST  FOURIER  TRANSFORM  PROCEDURES 


The  details  given  in  this  Appendix  apply  specifically  to  the  Theodorsen 
Garrick  mapping;  with  obvious  changes  in  notation,  they  also  are  valid  for  the 
derivative  form  of  the  transformation. 


Equations  2-20  and  2-21  contain  2N*-2  constants,  which  are  evaluated 
as  follows:  first,  they  are  satisfied  at  a  discrete  number  of  points,  denoted 
by<^^  : 


0K  * 


art  (.<-  1) 


K 


, 


CA-1) 


where  N  was  chosen  to  be  64,  in  the  present  case.  Thus: 

N-  1 

'Tc  •  *  Z  (^4  ^  i  ~  -i 

f  '  ' 


CA-2) 

(A- 3) 


Each  right-hand  side  now  contains  2N  coefficients.  The  correspondence  between 
either  of  these  and  the  Fast  Fourier  Transform  (FFT)  as  presented  in  Ref.  12 
is  given  by  the  next  two  equations:  consider  the  expression 


y  (1) 


2N-I 

L 

it 


i 


0,1  ;  W. 


ZN 


ZTT 

(A-4) 


where  values  Y(-)  are  real,  and  the  2/V  values  of  C(’l  are  in  general  complex, 
but  must  satisfy  the  following  redundancy  condition,  in  order  that  the 
values  be  real: 

CCn)  =  C(ZN-n)  ,  n  -  7,2,...,  /V  - ;  (A-S) 


C  (0)  and  CfA/)pure  real 

where  the  tilde  denotes  the  complex  conjugate.  When  these  conditions  are 
met,  Eq.  A-4  can  be  written  as 

^  (J.)  =  +  (“  f )  Cff(N)  +  2  c««  — in)  ^ 

■1“  o,  2A/-;  Ca-6) 

This  form  can  now  be  used,  in  conjunction  with  Eq.  2-20  or  2-21,  to  facilitate 
application  of  the  FFT  to  the  complex  form  given  in  Eq.  A-4. 


In  the  case  of  Eq.  A- 2,  values  of  the  coefficients  fl,  through 
and  6, through  are  found,  from  given  values  of ^  .  This  is  procedure 
4  of  Ref.  6,  which  takes  the  following  steps: 


1.  Set  X/tI)  -  Y(2^)  - 

X,(4)  =  Y(2-fcfO  - 


-fc  -  0,7,...,  A/-7 


2.  Set 

3.  Calculate  the  W  -point  Discrete  Fourier  Transform  of 

,  -ni 

/Ren)  -  +  c/R^in)  -  —  ^  Xc^)  ;  n  •  o,i, 


(A- 7) 


7  (A-8) 


/V-J 


(A- 9) 


4.  By  periodicity,  set  /R  (N)  “  /R  (o) 

5.  Apply  Eq.  34  of  Ref.  12,  in  order  to  extract  /fl,^«)and  from 

/R^(n)  .  .i  ^ 

.  ,  ..  *  n  •  0,  1 ,  y  CA-10) 

/R^{.n)  »  -I"  (N-i) 

Note  that  these  expressions  use  #)(rt>for  *1*0,  to  give  and 

for  ■  O,  7,...^  A//2 

6.  These  values  of  and  iR^(n)  then  give  C/n)  for  the  same  range: 

C(n.)  -  •“  *  0,1 ,  .  ,  N/Z  (A-11) 

N 

For  the  range  of  n  from  ■5-  7- 1  to  N-1  ,  use  Eq.  36  of  Ref.  12,  with  n  replaced 

*  ~n 

by  AJ- <7.  (and  noting  that  *  -  7  ): 

C(/a}  -  C  (2N-n) 


(A- 12) 


C,{A.)  s  ^  Cz  -  ft) 


Replace  rt  by  N  +  ia  : 


C(N-t-n)  a  cCzN-N-n)  *  C(N-rt') 


Thus  Eqs.  40  and  41  of  Ref.  12  are 

.  C(n.)  ^  c(N~rt.) 

r  ^ 

These  are  all  the  values  needed  for  and  . 


^  L  N/z 


(A- 19) 


3.  Find  from  Eqs.  42  and  43  of  Ref.  12: 


fi(N’-rt)  a  (/t)  •#  t /?^(n) 


n  ~  o,  A//2 


This  gives,  on  the  left-hand  side,  all  values  from  to  n*-N. 


4.  Calculate 


(A- 20) 


X(^)  '  Z  ^  ^  ^  ^  0,  N-1 


(A- 21) 


5.  Finally: 

-®o  -  ^ 

*  -k  0,1^  N-1 

-  S.  '  J-  [><(•*>]  J 

The  first  of  these  equations,  with  "^-o,  is  used  to  find  8*  . 


In  order  to  apply  these  formulas,  it  is  necessary  to  have  a  relation 


between  r  and  ^  : 


«  r(0^')  ,  K-1,Z,...,ZN  CA-23) 


which  is  foxmd  from  a  spline  fit  to  the  blade-surface  image  in  the  u?-  plane. 


The  discrete  Fourier  transform  and  its  inverse  are  given  by 


iR  (n)  - 

j_ 

N 

L  , 

/)  -  0^1 ,2,  N-1 

(A- 24) 

r" 

N- 

1  • 

X 

n 

L 

■i  =  0,  1 , 2,  ,  N-1 

(A- 25) 

n  « 

o 

0 

The  IMSL  routine  FFT2  evaluates  the  second  of  these,  i.e.,  it  returns 
X(^)  ,  given  the  values  /fl  (.n)  .  To  evaliiate  the  first  of  these,  FFT2  is 
used  with  input  /  N  ,  and  with  output  interpreted  to  bs  fi(n)  i 

this  is  Procedure  1  of  Reference  6: 


^  A/-I  ^  .  ni 

NAin)  -  2  X(^) 

i-o 


(A- 26) 


In  the  FORTRAN  version  of  these  and  other  procedures,  it  is  conve¬ 
nient  to  use  indices  that  begin  at  one,  rather  than  zero,  by  setting  ^ 

(FORTRAN  symbol  JP).  The  corresponding  table,  for  example  of  the  coefficient 


In  addition,  care  must  be  taken  with  the  argument /V-n. ;  for  example,  Eqs. 
are  written  as 


2ftl(NP)  =  ZCC(NP)  +  zee  (N  -hZ  -  ) 

Zf)Z(NP)  »  [ZCC(NP)  -  zee  (N*Z-  NP)'\ 


NP=  r,2, 


(A- 27) 


It  can  be  verified  that  the  quantity  N  +  2  -  NP  in  the  argtunents  above  pre¬ 
serves  the  correct  ordering  -  for  example  when  N  =  64,  and  n  =  0,  e( N~n) 
This  would  be  stored  in  ZCC(64  +  2  -  1)  =  ZCCf65). 


>  - 
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Appendix  B 

COMPUTER  PROGRAM  LISTING 


“RQGRAM  RAOESC INPUT, OUTPUT, GRID, TAPE5= INPUT r  ThPES  =  OUTPUT, 

*  TAPE7=GRID) 

PROGRAM  RAMES:  THE  lUES-LIUTERMOZA  CONFORMAL  TRANSFORMATION  FOR 
TURBOMACHINERY  CASCADES  (AIAA  JOURNAL, VOL.  5,1977,  PP  547  -  552) 
DOCUMENTATION  IS  GIUEN  INI 

N,  J.  RAE,  A  COMPUTER  PROGRAM  FOR  THE  IVES  TRANSFORMATION 
IN  TURBOMACHINERY  CASCADES,  CALSPAN  CORPORATION  REPORT  S275-A-3 , 
NOVEMBER  1980 

W.J.  RAE,  MODIFICATIONS  OF  THE  IVES  -  LIUTERMOZA  CONFORMAL- 
MAPPING  PROCEDURE  FOR  TURSOMACHINERY  CASCADES,  ASME  PAPER 
83-GT-llG,  MARCH  1983 

W.J.  RAE,  REVISED  COMPUTER  PROGRAM  FOR  EVALUATING  THE  IVES 
TRANSFORMATION  IN  TUR30MACHINERY  CASCADES,  CALSPAN  CORPORA~!ON 
REPORT  7177-A-l,  JULY  1983 

W.J.  RAE  AND  P.V.MARRONE,  RESEARCH  ON  TURBINE  FLCNFIELD 
ANALYSIS  METHODS,  CALSPAN  CORPORATION  REPORT  7177-A-3, 

JANUARY  1985 

IMPLICI-!-  REAL(A-H,a-Y)  ,COMPLEK(Z) 

LOGICAL  PNCHZA 

CaMMON/TGINTG/N,NPl , NP2 , N2 , NB2 , NB2P 1 , IP, IPMX, ITP, IWK , IMX,K JMX 
COMMON/TGCMPX/Zl , ZW2N , ZI , ZNN , ZA , ZCC , ZA 1 , ZA2 

C  ;MMON/TGDBLE/PBN  ,  OM  ,  CMM ,  ANGERR ,  A  ,  B ,  E ,  F  ,  THT  ,  :=H  I  ,  X  ,  Y  ,  BETA  ,  S  -  DR  ,  D I 
0  JMENSION  ZS(80)  ,ZP(80)  ,ZaMS(80)  ,Z0MP(80.> 

DIMENSION  RDS(aO) ,RDP(80) ,THS(aO) ,THP(BO) 

DIMENSION  RDSX(aO) ,RDPX(80) ,THSX(SO) ,THPX(BO) 

D: MENS  I ON  X( ISO) , Y( ISO) ,E( ISO) ,F( ISO) ,THT( ISO) , 

*  BETA< ISO) ,S( ISO) 

D- MENS I ON  PHI ( 130) , A(G5) ,B(G5) ,ZCC(S5) , DR ( 55 ) , D I ( G5 ) , 

*  ZA( 155) ,ZA1 ( 165) ,ZA2( 155) , 

•»  ZEE(  155)  ,ZFF(  1G5)  ,ZSaMGA(  165)  , 

*  ZX I ( 1 0 , 40 ) , ZEETA (10,40), ZETA (10,40), 2DMA (*0,40), 

ZFNL(  10,40)  ,ZXY(  10,40)  ,ZVEL(  10,40)  ,ZFNP(  10,40)  , 

*  I WX ( 7 ) 

DIMENSION  ITP( 100) , ID(3G) 

DIMENSION  XX (50, 20) , YY(50,20) 

DC  10  I  =  1,7 

10  IWK (I)  =  0 

NAMEL I  ST / INPUTS/  ANGERR , EX , G , H , I GOT , I LE , I TE , K MX , AL P , 

*  IMX,  LMX,  CM,  PNCHZA,  RTOL,  ZC,  ZLE,  ZN,  ZT ,  ZTZ, 
IPMX,  ISHEAR,K  JS,KJP,  :CIRC,NS 

-  SET  NAMELIST  DEFAULT  VALUES 


.'-V-l-.-; 


•  *  *’  V  H 


•  ’.■*  '/  V  ' 


L 

;CIRC=: 
ANGERR=0.024 
ALP  =  0. 
IG0T=0 


KMX  =  40 
LMX  =  10 
KJS  =  20 
KJP  =  20 
I  MX  =  400 
I SHEAR  =  1 
IPMX  =  1 
NS  =  40 
PNCHZA=. FALSE. 

RTOL  =  3. 

OM  =  1  . 

C 

C  note:  EXr  G,  H,  ZCr  ZLEr  ZN ,  ZT,  ZTE  HAOE  ND  DE-AL'LT  MALUES 

C  (ZC  IS  NOT  NEEDED  IF'  I60T  =  0) 

C 

READ(5, 103) ( ID( I) r 1=1 r3G) 

103  FCRMAK  iaA4) 

r' 

PI  =  4,0  *ATAN(1.0  ) 

TPI  =  2.0  >PI 
P32  =  PI/2. 

ZPI  =  CMPLX(PI,0.0  ) 

ZT  =  CMPLXd.O  ,0.0) 

ZERO  =  CMPLX(0.0  ,0.0) 

.“MGA  =  CMPLX  (  + 1 . 0  ,  0 . 0  ) 

ZMGB  =  CMPLX (-1.0  ,0.0) 

Z'- PS  =  CMPLX  ( 0.  ,  1  .E-08) 

Z-'.G  =  CMPLXd  .ElO,  1  .ElO) 

C  READ  PNCHZA  =  T  IF  ( ZA  (  L )  ,  L=  1  ,  LMX )  IS  TO  BE  PUNCHED  FOR  ALL  '.'SLUES 
C  OF  K,  OTHERWISE  READ  PNCHZA=F 

C  I  GOT  =  1  IF  THE  UALUE  OF  ZC  IS  KNOWN.  OTHERWISE,  I  GOT  =  0. 

C  KMX  AND  LMX  ARE  THE  GRID  SIZES  IN  THE  FINAL  TRANSFORMED  PLANE. 

C  KJS  AND  KJP  ARE  THE  NUMBERS  OF  POINTS  ON  THE  SUCTION  AND  PRESSURE 
C  SIDES  AT  WHICH  PAIRS  OF  BLADE  COORDINATES  WILL  BE  INPUT. 

C  ILE  =  0  OR  1  FCR  A  SHARP  OR  ROUNDED  LEADING  EDGE,  RESPECT  I '.'ELY 
C  ITE  =  0  OR  1  FOR  A  SHARP  OR  ROUNDED  TRAILING  EDGE,  RES  PECT I '.'SLY . 

C  FOR  A  SHARP  LEADING  EDGE  <ILE=0)  ZN  MUS^  EQUAL  ZLE 
C  FOR  A  SHARP  TRAILING  EDGE  (ITE=0)  ZT  MUST  EQUAL  ZTE 

C  IMX  IS  THE  MAXIMUM  NUMBER  OF  ITERATIONS  ALLOWED  FOR  THE  PH  I /■'’HE'^A  ( /S ) 
C  ITERATIONS 

C  IPMX.NE..  CAN  BE  USED  TO  DISPLAY  THE  '.'ALUE5  OF  THETA/S  DURING  THE 
C  PH:/THETA(/S)  iterations,  at  the  iteration  numbers  READ  INTO  THE 
C  I  TP  ARRAY  BELOW. 

C  ISHEAR  =  0  GI'.'ES  AN  ORTHOGONAL  GRID.  THE  LINES  K  =  1  AND  K  =  KMX,  WHICH 

C  START  AT  THE  TRAILING  EDGE,  DO  NOT  GO  TO  DOWNSTREAM  INFINITY. 

C  ISHEAR  =  1  SHEARS  THE  GRID  UNIFORMLY!  I-N  THIS  CASE,  THE  K  =  1  AND  KMX 

C  LINES  DO  GO  TO  DOWNSTREAM  INFINITY. 


B-2 


.'•V 


OM  IS  A  RELAXATION  FACTOR  USED  IN  THE  PHI/THETA  MAPPING: 

-  FOR  I CIRC  =  0,  USE  0.1,  OR  A  SMALLER  UALUE  IF  THE  A  AND 
ITERATIONS  FAIL  TO  CONVERGE 

-  FOR  ICIRC  =  1,  DM  =  1.0  SHOULD  BE  SATISFACTORY. 

ANGERR  IS  THE  ANGULARdN  RAD  I ANS )  ( /ARC  LENGTH)  TOLERANCE 
FOR  THE  PHI/THETA(/S)  TRANSFORMATION. 

A  REASONABLE  VALUE  IS  .01 

RTOL  IS  THE  TOLERANCE  FOR  THE  MAX/MIN  RADIUS  RATIO  IN  "HE 
OMEGA  PLANE 

ALP  IS  THE  ANGLE  OF  INCIDENCE,  IN  DEGREES  (NEGATIVE),  IN 
THE  ETA-PLANE. 

ICIRC  =  0  FOR  THE  THEODORSEN-GARR I CK  TECHNIQUE 
ICIRC  =  1  FOR  THE  BAUER  ET  AL  TECHNIQUE 

NS  IS  THE  NUMBER  OF  STEPS  TO  BE  USED  IN  THE  SIMPSON 'S-RULE 
INTEGRATION  OF  D(SMALL  OMEGA ) /D ( ZETA )  WHEN  ICIRC  =  1. 

-  READ  NAMELIST  INPUT  DATA 

READ (5, INPUTS) 

ITP(i)  =  999 

IF( IPMX.NE. 1 )  READ(5, 1 02 ) ( ITP ( I P ) , IP=1 , IPMX) 

102  FORMAT (2014) 

IP  =  1 


KMXH= (KMX+1 ) /2 
KMXL=KMXH 

IF(M0D(XMX,2) .NE.O) 
lOE  =  M0D(KMX,2) 
KOUNT  =  0 


.  =  KMXH-1 


DO  310  K=l,50 
Dr'  311  L=1,20 
XX ( K , L )  =  0 . 
YY<K,L)  =  0. 
CONTINUE 
CONTINUE 


INPUT  VARIABLES  EX,  G,  H,  ZLE,  ZN,  ZT,  ZTE, 

THE  FORTRAN  STATEMENTS  IN  SUBROUTINE  SHAPE, 

AND  THE  VARIABLES  LISTED  IN  COMMON  BLOCK  GEOM  (IF 
ARE  ALL  SPECIFIC  TO  THE  BLADE  SHAPE  BEING  USED. 


IS  BEING  USED) 


CALCULATION  OF  THE  BLADE  SHAPE 


D  =  CABS(ZTE-ZLE) 

CALL  SHAPE ( D , H , G , EX , ZP , ZS , K JS , K JP ) 
SO  =  SQRT<H-»-H+G*G) 

ZDA  =  CMPLX(G/SG,-H/SG) 
ZGAMMA=CONJG(ZDA) 

A:_P  =  ALP*TPI/3B0. 

ZALP  =  CMPLX(COS(ALP) ,SIN( ALP) ) 

ZAM  =  CDNJG(ZALP) 

ZAP  =  ZALP 
k'RITE(S,207) 

207  FORMATdHl; 


13X,  'N',13>(,  'X',13X, 


|.;RITE(B,20B)  (  ID(  I )  F  1  =  1 ,3B> 

20B  FORMAT ( 30X, 18A4) 

KRITE(S,240)  DfHfGfEXfSG 

240  FORMAT (//I OX F  'BLADE-GEOMETRY  PARAMETERS  AREI'f 
//5X,  'ABS(ZTE-ZLE)  ='f 

*  fio.Sf'  h  =  'fFio.Sf'  g  =  'fFio.Sf'  ex  = 

•»  FIO.Sf  '  SLANT  GAP  ='fF10.5f//) 

WRITECSf INPUTS) 

KJLE  =  KJP  +  2 
KJMX  =  KJLE  +  KJS  +  1 
WRITE (Bf 203) 

203  FORMAT (  9Xf  'BLADE  COORDINATES t' f 

*  //3Xf 'SUCTION  SIDE' f/3Xf 'KJ'f7Xf 'S' 

ZZXY  =  ZGAMMA*2LE 
U;RITE(Bf270)  KJLE  fZLEf  ZZXY 

270  FORMATf I5f 1P4E14.5) 

DO  S4  K  =  IfKJS 
KJ  =  KJLE  +  K 
ZZXY  =  ZGAMMA*ZS(K) 

KRITE(Bf270)  KJfZS(K) fZZXY 

64  CONTINUE 

ZZXY  =  ZGAMMA*ZTE 
WRITE (Bf270)  KJMXfZTE  ,ZZXY 
W;-' I  TECS,  271  ) 

271  F0RMAT(//8Xf 'PRESSURE  SIDE ' , /3Xf 'K J ' , 7Kf 'S ' , IGX, 'N ' , 

*  13Xf  'X'f13Xf  'Y',/) 

ZZXY  =  ZGAMMA*ZLE 
WR I TECS, .270)  K JLE , ZLE , ZZXY 
DC'  SI  K  =  1  F  K  J P 
KJ  =  KJLE  -  K 
ZZXY  =  ZGAMMA*ZPCK) 

WRITE CSf 270)  K J , ZP C K ) , ZZXY 

SI  CONTINUE 
KJ  =  1 

ZZXY  =  ZGAMMA*ZTE 
WRITECSfZ70)  KJfZTEfZZXY 


ZDN  =  CMPLXCHfG) 

ZZ  =  C2T-ZN)/ZDN 

CMI  =  PI*EX*CG*REALCZT-ZN)-H*AIMAGCZT-2N) )/CSG»SG) 
XA  =  PI#EX*CH*REALCZT-ZN)+G*AIMAGCZT-ZN) )/ CSG*SG) 

R  =  EXPC-CHI) 

ZPLUS  =  CMPLXCR*COSCXA) f-R*SINCXA) ) 

R  =  l.O/R 

ZMINUS=  CMPLXCR*COSCXA) fR*SINCXA) ) 

DO  20  K  =  1 fKJS 
Z.'T  =  CZSCK)-ZT)/ZDN 
ZETAIS  =  2PI*2ZT 
2ETA2S  =  CSINCZETAIS) 

ZZN  =  CZSCK )-ZN>/ZDN 
ZETASS  =  ZPI-^ZZN 
ZeTA4S  =  CSINCZETA3S> 

ZFS  =  ZETA2S/ZETA4S 


1^' WWmyf  JUikUL'A'  H^iLlLUIL" 


ftA.  v:\yv  AT 


p 

t? 

p 


c 

c 

c 

c 


c 

c 

c 

c 


RDS(K)  =  CABS(ZFS) 

THS(K)*  ATAN2(AIMAG(2FS) ,REAL(ZFS) ) 

20  CONTINUE 

DO  24  K  =  1,KJP 
ZZT=  (ZP(K)-2T)/ZDN 
ZFTAIP  =  ZPI*ZZT 
ZFTA2?  =  CSIN(2ETA1P) 

ZZN  =  (ZP(K >-ZN) /ZDN 
ZETA3P  =  ZPI»2ZN 
ZETA4P  =  CSIN(2ETA3P) 

ZFP  =  ZETA2P/ZETA4P 
RDP(K)  =  CABS (ZFP) 

THP(K)=  ATAN2(AIMAG(ZFP) .REAL (ZFP) > 

24  CONTINUE 

NOW  ADD  THE  LEADING-  AND  TRAILING-EDGE  POINTS.  AND  STORE  THE  G(Z) 

ARRAY  AS  E(KJ)*EXP(I*THT(KJ) ) .WHERE  KJ=1,KJMX  AS  YOU  GO  FROM  TE  AROUND 
THE  PRESSURE  SIDE  TO  THE  LE  (KJ=KJLE)  AND  THEN  ALONG  THE  SUCTION  SIDE 
SACK  TO  THE  TE  AGAIN  (KJ=KJMX). 

IF( ITE.EQ. 1 )  GO  TO  13 
E( 1 )=0.0 
THT( 1 ) =0.0 
GO  TO  14 

13  ZETA2=CSIN(ZPI*(ZTE-ZT)/2DN) 

ZFTA4=CSIN(2PI*(2TE-ZN)/2DN) 

2FP=ZETA2/ZETA4 

E( 1 )=CA3S(ZFP) 

THT( 1 )=ATAN2(AIMAG(ZFP) ,REAL(ZFP) ) 

14  IFdLE.EQ.l)  GO  TO  15 
E(KJLE)  =  l.OE+OS 

THT(KJLE)  =  0.5*(THP( 1 )+THS( 1 ) ) 

GO  TO  IG 

15  ZETA2=CSIN(2PI*(2LE-ZT) /ZDN) 

ZETA4=CSIN(2PI*(ZLE-ZN)/ZDN) 

ZHP=ZETA2/ZETA4 

E(KJLE)  =  CABS(ZFP) 

THT(KJLE)  =  ATAN2(AIMAG(ZFP) .REAL(ZFP) ) 

IS  DO  17  K  =  l.KJP 
KJ  =  KJLE  -  K 
E(KJ)=RDP(K) 

THT(KJ)=THP(K) 

17  CONTINUE 

DO  18  K  =  l.KJS 
KJ  =  KJLE  +  K 
E(KJ)=RDS(K ) 

IB  THT(KJ)=THS(K) 

E(KJMX)  =  E( 1 ) 

THT(KJMX)  =  THT( 1 ) 

NOW  ADJUST  THE  BRANCHES  OF  G(Z)  SO  AS  TO  BE  CONTINUOUS  ACROSS 
THE  CUT  (ALONG  THE  NEGATIVE  REAL  AXIS)  OF  THE  ATAN2  FUNCTION. 

251  BR=0.0 
KA  =  2 

KB  =  KJLE  +  1 
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IF(ITE.EQ.O)  KA  =  3 
IF(THT(KA-1 ) .GT.O. >  BR  =  -1. 

Pa=THT(KA-l ) 

DO  321  KJ=KA,KJLE 
C.-;G  =  THT(KJ)-PO 
pa=THT(KJ) 

IF(ABS(CHG) .LE.PI )  GO  TO  321 
IF(CHG.GT.?I )  8R=BR-1.0 
IF(CHG.LT.-PI)  BR=BR+1.0 
THT(KJ)=THT(KJ)+BR*TPI 
( ILE.EQ. 1 )  GO  TO  32G 
ER  =  BR  +  1 . 

PQ  =  THT(KB> 

THT(KB)  =  THT<K3)  +  BR*TPI 

KB  =  KB  +  1 

DO  327  KJ  =  KB,KJMX 

CHG  =  THT(KJ)-PO 

PO  =  THT(KJ) 

IF(ABS(CHG) .LE.PI )  GO  TO  327 
IF(CHG.GT.PI)  BR  =  BR  -  1. 

IF(CHG.LT.-PI )  3R  =  BR  +  1. 

THT(KJ)  =  THT(KJ)  +  BR^TPI 
DO  323  K=1,KJP 
THP(K)  =  THT(KJLE-K) 

CONTINUE 

Du  322  K  =  i.KJS 
THS(K)  =  THT<KJLE+K) 

CONTINUE 

DO  25  K  =  1,KJS 
ARS  =  £X*THS(K) 

RS  =  RDS(K)**EX 
RDSX(K)  =  RS 
THSX(K)  =  ARS 
CONTINUE 
DO  2B  K  =  1 ,KJP 
ARP  =  EX*THP(K) 

RP  =  RDP(K)**EX 
RDPX<K)  =  RP 
THPX(K)  =  ARP 
CONTINUE 
WRITE (B, 241 ) 

FORNAT( //I  OX,  'BLADE-SURFACE  IMAGES  IN  THE  G  -  PLANE 
<•  '  SINES)  AND  CAP  OMEGA  PLANE  ( G**l /KAPPA  >  ,AND', 

^  //  9X,'  RADII  AND  ANGLES  USED  IN  SELECTING  THE  P 

f  '  BRANCHES  OF  THE  RATIO  OF  SINE  FUNCTIONS  ARE!', 
^  //JX,  KJ  '  , 15X,  'G' ,25X,  'CAP  OMEGA ', 13X ,' R ', 1 1 X ,' T 

^  3X,  'R**EX  '  ,7X,  'THETA*EX  ' , ) 

XS=E(  l)»COS(THT(  1)) 

YS=E(  1)*SIN(THT(  1)) 

RDTX  =  E(1)**EX 
THTX  =  EX*THT(1) 

US  =  RDTX*COS(THTX) 

OS  =  RDTX*SIN(THTX) 


(RATIO 


PROPER  ■'  , 


THETA  '  r 


: 


.V.*. 

."I-/  -y*. 

.  u”,  .>  . 
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-  V  '.\N 


k:?  I TE  (  S  ,  325 )  XS  r  YS  ,  US ,  VS  ,  E  ( 1  )  ,  THf  ( 1  )  ,  f?DTX ,  THTX 
232  FORMAK  I5r  lPaE14.5) 

324FC'RMA7('  LE  ',iPSEi4.5) 

325  FORMAK'  TE  ',1P8E14.5) 

KJLM  =  KJLE  -  1 
DO  B5  KJ  =  2rKJLM 
K  =  KJLE  -  KJ 
X=  =  RDP<K)*COS(THP(K  )  ) 

YP  =  RDP(K)*SIN<THP(K) ) 

UP  =  RDPX(K)*COS<THPX(K) ) 

=  RDPX(K)*SIN(THPX(K) ) 

XRITE(S,232)KJ,XP,YPfUP,VP,RDP(K) ,THP(K) ,RDPX(K) ,THPX(K; 

St-  CONTINUE 

XS  =  E(K JLE)*CaS(THKKJLE)  ) 

Y3  =  E(KJLE)-»SIN(THT(KJLE)  ) 

RDLX  =  E(KJLE)**£X 
THLX  =  EX*THT(KJLE) 

US  =  RDLX*COS(THLX) 

VS  =  RDLX*SIN(THLX) 

WR I TE ( B , 324 ) XS , YS  r  US , VS , E ( K  JLE ) r  THT ( K  JLE ) r  RDLX , THLX 

KJMXM  =  KJMX  -  1 

KJLP  =  KJLE  1 

DO  31  KJ  =  KJLP, KJMXM 

K  =  KJ  -  KJLE 

XH  =  RDS<K )*COS(THS(K  )  ) 

YS  =  RDS(K )*SIN(THS(K ) ) 

US  =  RDSX(K)-»COS<THSX(X)  ) 

VS  =  RDSX(K )»SIN(THSX(K )  ) 

kR I TE ( S  r  232 )  K  J , XS , YS , US , VS , RDS  <  K ) , THS ( K ) , RDSX ( X ) , th5X ( K ) 

51  CONTINUE 

X  P  =  E ( K  J  MX ) *COS ( THT ( K  J  MX ) ) 

YP  =  E(KJMX)»SIN<THT(KJMX) ) 

RDTX  a  E(KJMX)**EX 
THTX  a  £X*THT(KJMX) 

UP  a  RDTX*CaS<THTX) 

V?  a  RDTX*SIN(THTX) 

WRITE(BF325)XPrYP,UP, VP,E(KJMX) , THT (KJMX) ,RDTX,THTX 
WPITE(Br242)ZPLUS,ZMINUS 

242  FaRMAT(//10X, 'POINTS  AT  INFINITY  ARE  LOCATED  IN  THE  CAP  OMEGA', 

♦  '  PLANE  at:  ' , 

>  //15X,  'PLUS:  ' , 1P2E15.4,  '  MINUS :  ' , 2E 1 5 . 4 , / / > 

C 

C  .DETERMINATION  OF  ZC  SUCH  AS  TO  MINIMIZE  THE  RATIO  RMAX/RMIN  IN  THE 
C  L.C.  OMEGA  PLANE 

M=1 

ZE  a  (ZMCiA-ZMCiB)  /  (ZPLUS-ZMINUS) 

Z!-  a  (ZMGA#ZMINUS-ZMGB*2PLUS)/(ZPLUS-ZMINUS) 

2G  a  (ZMGA^^ZPLUS-ZMG3^^ZMINUS)/(ZPLUS-ZMINUS) 

WRITE (G, SOI ) 

GOl  FORMAT(  'i!TER',llX,  'ZD',22X, 'Z8',Z2X,  'ZC',iaX,  'ZOMSTR ' , 19X ,  'ZNTRD' 
1 /i3X,  'RMIN  '  ,8X,  ’RMAX'fTX,  'RATIO'/'  (ZAiKJ) ,KJ  =  1,KJMX)  ') 

250  ITERal 

RAT 1 0  =  0.0 

IF( IGOt.EQ. 1 )  GO  TO  GO 
C 

C  FDR  A  FIRST  GUESS,  USE  ZC= ( -1 . 0 , +1 . 0 ) 

C 


ZC=CMPLX(-1 .0, i .0) 

GO  ZB  =  (ZMGA*ZPLUS-ZMGB*ZMINUS-ZC-»(ZMGA-ZMGB)  )  /  (ZPLUS-ZMI 'iUS) 
ZD  =  ( ZMGB»ZPLUS-ZMGA*ZM  INUS+  ( ZrGA-ZMGB  )  /ZC  >  /  (  ZPLUS-ZM I  NL'S  ; 
GB  CONTINUE 

DO  75  K  =  1,KJS 
RS  =  RDSX(K) 

ARS  =  THSX(K) 

ZOMS(K)  =  CMPLX(RS*CQS(ARS) .RS*SIN(ARS) ) 

ZHMSiK)  =  (ZD-ZB»Z0MS(K)/ZC)/<Z1-Z0MS(K)/ZC) 

75  CONTINUE 

DO  85  K  =  1,KJP 
RP  =  RDPX(K) 

ARP  =  THPX(K) 

ZOMP(K)  =  CMPLX(RP*COS(ARP) ,Ro*SIN(ARP) ) 

ZOMP(K)  =  (ZD-ZB*Z0MP(K)/ZC)/(Z1-Z0MP(K5/ZC) 

85  CONTINUE 

I?^(  ILE.EQ.  1  )  GO  TO  11 
ZOMLE  =  ZB 
GO  TO  12 

11  RD  =  E(KJLE) 

TH  =  TH'^(KJLE) 

RS  -  RD+'ff’EX 

TK  =  EX*TH 

ZOMLE  =  CMPLX(RS*COS(TH) ,RS*SIN(TH) ) 

ZOMLE  =  (2D-2B*ZCMLE/ZC>/(21-Z0MLE/ZC) 

12  IF( ITE.EQ. 1 )  GO  TO  32 
ZOMTE  =  ZD 

ZA(i)  =  ZOMTE 
GO  TO  33 

32  RD=E(1) 

Th=THT( 1 ) 

RS=RD**EX 

TH=EX*TH 

ZOMTE=CMPLX(RS*COS(TH) ,RS*SIN(TH) ) 

20MTE=(ZD-ZB*Z0MTE/ZC)/ (Zl-ZOMTE/ZC) 

ZA( 1 ) =ZOMTE 

33  DO  7G  KJ  =  2,KJLM 
K  =  KJLE  -  KJ 

7G  ZA(KJ)  =  ZOMP(K) 

ZA(KJLE)  =  ZOMLE 
DO  77  K  =  1 ,KJS 
KJ  =  KJLE  +  K 
77  ZA(KJ)  =  ZOMS(K) 

ZA(K JMX)  =  ZA(  1  ? 

ZNTRD  =  CMPLX(O.OrO.O) 

AREA  =  0.0 
RMIN=CABS(ZA( 1 ) ) 

RMAX=RMIN 
ZMAX=ZA( 1 ) 

2MIN=ZA( 1 ) 

KJMXX  =  .1 
KJMN=1 

DO  78  KJ  =  2,KJMX 

DAREA  =  ABS(REAL(ZA(KJ-1 ) )*AIMAG(ZA(KJ) )-REAL(ZA(KJ) 

*  AIMAG(ZA(K J-1 ) > ) /2.0 
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c 

C  SE1 
C 


7!3R  =  (ZA(KJ-1)+ZA(KJ)  )/3.0 
ZNTRD  =  ZNTRD  +  ZBR*DAREA 
Ri»B5  =  CABS(ZA(KJ)  ) 

I'-  (RABS.GE.RMIN)  GO  TO  79 
R,'1IN  =  RABS 
Z'-IIN  =  ZA(KJ> 

K JMN=K J 
GO  TO  78 

IF(RABS.LE.RMAX)  GO  TO  78 

Rr,AX  =  RABS 

ZHAX=ZA(KJ) 

KJMXX  =  KJ 

AREA  =  AREA  +  DAREA 

RATIQ=RMAX/RMIN 

ZNTRD  =  ZNTRD/AREA 

ZOMSTRs  ZC* ( ZNTRD-ZD ) / ( ZNTRD-ZB ) 

K!RITE(GrG02)  ITER  ,  ZD  ,  ZB  ,  ZC  ,  ZOMSTR , ZNTRD  ,  RM I N  ,  RMAX  ,  RATI 0  , 

L  (2A(KJ) ,K J=1 ,KJMX) 

FORMAT (/ 15 r lP10E12.4/Ei7.4,2E12.4/(10E13.5) ) 

I-(RATIO.LT.RTOL)  GO  TO  G3 
I- ( IGOT.EQ. 1 )  GO  TO  G3 
r'ER  =  ITER  +  1 
Ii='(  ITER.LE.30)  GO  TO  S2 
:-:RITE(5,204) 

FaRMAT(///10X,  'TOLERANCE  SPECIFIED  FOR  RMAK/RMIN  NOT  MET  IN'r 
»  '  30  ITERATIONS') 

S  ' OP 

CONTINUE 

IF(M.Ea.2)  GO  TO  SB 
K  =  2 

2DS=1 . 1*ZMIN 
K J=K JMN 
RD  =  E(i<  J  )»*EX 
TH=EX*THT(K J ) 

ZOM=CMPLX(RD*CaS(TH) ,RD*SIN(TH) ) 

ZC=  ( ZOM*  ( ZDS-ZG )  +ZE )  /  ( ZDS+ZF-ZE-^ZOM ) 

GO  TO  GO 
M=1 

2DS=0.9  *ZMAX 
KJ=KJMXX 
GO  TO  G7 
IGOT  =  1 

WR I TE ( B  7  208 )  ZD , ZB , ZC  »  ZNTRD , ZOMSTR 
FORMAT <//10X 7  'CONSTANTS  FOR  MAPPING  FROM  '7 
»  'CAP  OMEGA  -  PLANE  TO  SMALL  OMEGA  -  PLANE  ARE  '  7 

//20X7'A  =  '  7  IP2E2O.57/2OX7  '3  =  '  7  1P2E20. 57/20X7  'C  = 

»•  iP2E20,5  7//20X7  'ZNTRD  =  '  7 1  P2E20 . 5  7 /20X  7  'ZOMSTR  =  ',:iPZE20.5) 

r  UP  THE  ARRAYS  OF  THETA  AND  LN(R) 


DO  41  KJ  =  i7KJMXM 

X(KJ)  =  ATAN2( AIMAG(ZA(K J ) ) 7REALCZA(K J >  )  ) 
Y(KJ)  =  CABS(ZA(KJ) ) 

X(KJMK)  =  X(l)  +  TPI 
Y(KJMX)  =  Y(l) 


c 

C  NOW  ADJUST  THE  ARGUMENTS  OF  THE  THETA  ARRAY,  SO  AS  TO  BE  CONTINUOUS 
C  ACROSS  THE  BRANCH  CUT  (ALONG  THE  NEGATIVE  REAL  AXIS)  OF  THE 
C  ATAN2  FUNCTION.  . THIS  ADJUSTMENT  ASSUMES  THAT  THE  CONTOUR  IS 
C  TRAVERSED  IN  A  COUNTERCLOCKWISE  DIRECTION. 

r 

8R  =  0.0 
PO  =  X'  ( 1) 

KOUNT  =  0 

DO  410  KJ  =  ZrKJMXM 
CHG  =  X(KJ)  -  PO 
IF(CHG.GT.O. )  GO  TO  409 
IF(CHG.LT.-PI )  GO  TO  408 
KC'JNT  =  KOUNT  +  1 
GO  TO  409 

408  BR  =  BR  +  1 . 

409  CONTINUE 
PO  =  X(KJ) 

X(KJ)  =  X(KJ)  +  BR*TPI 

410  CONTINUE 

IF(KOUNT.GT.O)  WRITE(G,407)  KOUNT 
407  F0RMAT(//5X,  'warning:  THERE  ARE', 13,'  NEGATIVE  INCREMENTS', 

*  '  IN  THE  ANGULAR  VALUES  OF  SMALL  OMEGA  ON  THE  3LADE  SURFACE', 
fr  //lOX,  'CHECK  WHETHER  THE  ZS  AND  ZP  ARRAYS  ARE  INTERCHANGED', 

//iOX,'OR  WHETHER  R  VS.  THETA  IS  MULT  I PLE-VALUED ,  IN  WHICH', 

'  case:  '  , //lOX,  '  -  IF  THE  OPTION  ICIRC  =  0  IS  BEING  USED,', 

*  '  THE  NEGATIVE  INCREMENTS  MUST  BE  REMOVED .',// 15X ,  'TRY  A', 

'  SMALLER  VALUE  OF  R-'-QL.', 

*  //lOX,'  -  IF  ICIRC  =  1,  NEGATIVE  INCREMENTS  ARE  ALLOWED:', 

■»  '  HOWEVER,  IF  THE  PHI/S  ITERATIONS  FAIL  TO  CONVERGE ',// i5X , 

•»  '  TRY  A  SMALLER  VALUE  OF  RTOL  OR  A  LARGER  VALUE  OF  ANGERR ' , 

*  '  OR  BOTH. ',//) 

C 

RMIN  =  10.0 

RMAX  =  0.0 

DO  49  K  =  1,KJMX 

IF( Y(K ) .LT.RMIN)  RMIN  =  Y(K) 

49  IF(Y(K ) .GT.RMAX)  RMAX  =  Y(K) 

WARSCH  =  SQRTIRMAX/RMIN)  -  1.0 
IF (WARSCH.lt. 0.3  )  OM  =  1.0 
WRITE(G,202) 

DO  43  KJ  =  1,KJMX 

ZSOMGA(KJ)  =  CMPLX(Y(KJ)-»COS(X(KJ)  )  ,Y(KJ)*SIN(X(KJ)  )  ) 

43  Y(KJ)  =  ALOG(Y(KJ) ) 

C 

C  2S0MGA  NOW  CONTAINS  KJMX  VALUES  OF  SMALL  OMEGA  ON  THE  BLADES 
C 

C  USE  FFT  TO  FIND  THE  "FOURIER  COEFFICIENTS  IN  THE  SMALL-OMEGA/ZETA 

c  mapping: 

C 

ZI  =  CMPLX(0.0, 1 .0) 

N  =  S4 


P8N  =  PI/FLOAT(N) 

M32  =  N/2 

NB2P1  =  NB2  +  1 

ZiMN  =  CMPLX  ( FLOAT  ( N  )  ,  0 . 0  ) 

ZW2N  =  CMPLX(COS(PBN) ,SIN(PBN) ) 
C’lM  =  1.0  -  DM 
DC  298  11=1,160 
BETAdI  )=0. 

S(ir)=0. 

298  CONTINUE 


FOR  ICIRC  =  1,  SET  UP  THE  CONSTANTS  FOR  THE  LN(R),  THETA  SPLIi 
IF(ICIRC.EQ.O)  CALL  C ISPLN ( Y , X , E , F , K JMX , 1 , 1 28 , 1) 


(ENTERED  WHEN  ICIRC=1)  COMPUTES  THE  BPTA  Awn  <5 
COORDINAiES,  AND  SET  UP  THE  CONSTANTS  FOR  THE  BETA~AND  S  SPL 


IF  (ICIRC. NE.l)  GO  TO  297 
ZAl ( 1 )  =  CMPLX(0. ,0. ) 

Dl  289  KJ  =  2,KJMX 

289  ZAl(KJ)  =  ZAKKJ-l)  +  ABS  (  ZSOMGA  ( K  J  ) -ZSOMGA  (  K  J- 1  )  > 
CALL  ZISPLN(ZS0MGA,ZA1 ,ZEE,ZFF,KJMX, 1 , 128,1 ) 

CALL  ZISPLN(ZS0MGA,2A1 ,2EE,ZFF,KJMX,4,12S, 1 > 

DO  290  KJ  =  1,KJMX 
BETA(KJ)  =  REAL(ZEE(KJ) ) 

290  S{KJ)  =  REAL(2FF(KJ) ) 

295  CALL  CISPLN(BETA,S,E,F,KJMX, 1 , 128,2) 

297  CONTINUE 


WRITE(6,243) 

243  F0RMAT(3X,  'BLADE-SURFACE  IMAGE  IN  THE  SMALL-OMEGA  PLANE:' 

*  /3X,  'KJ',BX,  'REAL'.lOX,  'IMAG',12X,  'R',9X,  'THETA', 

*  14X,  'S  '  , 14X,  'BETA ' , / ) 

DO  51  KJ  =  1,KJMX 

WRITE (6, 299)  KJ,ZA(KJ) ,EXP(Y(KJ) ),X(KJ),S(KJ) ,3ETA(KJ) 

299  FORMAT( 15, 1PGE14.5) 

5i  CONTINUE 


IF( ICIRC. EQ.O)  CALL  THDGRK 

IF( ICIRC. EQ. 1 )  CALL  BAUGRK ( RMAX , RMIN ) 


FIND  ZETAA  AND  ZETAB 

ZABST  =  ZMGA 
Z9BST  =  ZMGB 
RABST  =  1.0 
RBBST  =  1.0 
DRR  =  0.1 

DTH  =  10.0*PI/180.0 
THA  =  -9.0*DTH 
DO  21  I  =  1,9 
R  =  . 19  +  DRR*FLOAT( I-l ) 


V* 


D'l  22  J  =  1,19 

TH  =  DTH*FLOAT( J-1 )  +  THA 

ZTAGS  =  CMPLX(R*CDS(TH)  ,R*SI!M(TH)  ) 

IF  (ICIRC.EQ.O) 

CALL  CMETA(A.B,ZMG,ZTAGS.ZTANSR,S5,1.0E-00, 1  ) 

IF( ICIRC.EQ. 1 ) 

*r:ALL  OMDZETACZOMTE.Zl  ,  DR,  DI  ,ZMG,  ZTAGS,  ZTANSR,G5r  i  .OE-00, 1 ,2CO 
RA  =  CABS(ZMG-ZMGA) 

IF<RA.GT.RABST)  GO  TO  23 
RABST  =  RA 
ZABST  =  ZTAGS 

21  ZTAGS  =  -ZTAGS 
IF  (ICIRC.EQ.O) 

4  CALL  QMETA(A,B,ZMG,ZTAGS,ZTANSR,S5, 1 .OE-00, 1 ) 

IF( ICIRC.EQ. 1 ) 

♦CALL  0MD2ETA(ZQMTE,21 ,DR,DI ,ZMG,ZTAGS,ZTANSR,S5, 1 .OE-00, 1 ,20) 
R3  =  CABS(ZMG-ZMGB) 

I(^  (RB.GT.RE8ST)  GO  TO  22 
RBBST  =  RB 
ZBB5T  =  ZTAGS 

22  CONTINUE 
?1  CONTINUE 

ZTAGS  =  ZABST 
=  0 

IF  (ICIRC.EQ.O) 

♦  CALL  0META(A,B,ZMGA,ZTAGS,ZTANSR,B5, l.OE-OSrM) 

!*='(  ICIRC.EQ.  1  ) 

♦CALL  0MD2ETA(ZaMTE,Zl ,DR,DI , ZMGA , ZTAGS , ZTAN5R , B5 , 1 .OE-05,M,20> 
IF(M.NE.5)  GO  TO  2G0 
WR1TE(G,2G1)  ZTAGS, RABST 

2G1  F0RMAT(//5><,  'OMETA  FAILED  TO  CONVERGE  FOR  ZETA  A.’', 

*  /lOX, 'ZTAGS  =  ',lP2Ei3.5,'  RABST  =',E13.5) 

STOP 

230  CONTINUE 

ZGTAA  =  ZTANSR 
ZTAGS  =  ZBBST 
=  0 

IF  (ICIRC.EQ.O) 

♦  CALL  OMETA ( A, B,ZMGB, ZTAGS, ZTANSR, B5, 1 .0E-05,N) 

IF( ICIRC.EQ. 1 ) 

♦CALL  0MDZETA(Z0MTE,Z1 ,DR,DI , ZMGB , ZTAGS , ZTANSR , G5 , 1 .0E-05,M,20) 
IF(M,NE.5)  GO  TO  2G2 
l•:RITE(G,2G3)  ZTAGS, RBBST 

2G3  F0RMAT(//5X,  'OMETA  FAILED  TO  CONVERGE  FOR  ZETA  Bl  ' . 

*  /lOX,  'ZTAGS  =  ',1PZE13.5,'  RBBST  =',E13.5) 

5 'OP 

2G2  CONTINUE 

ZETAB  =  ZTANSR 

r* 

u 

C  FIND  GAMMA,  ALPHA,  BETA,  AND  S  FOR  MAPPING  TO  ETA  -  PLANE 
C 

A?  =  CABS(ZETAA  +  ZETAB) 

AM  =  CABSCZETAA  -  ZETAB) 

A3  =  CABS(ZETAA^ZETAB) 


.‘  V 


■.  .vV 


.  **-  ***  **• 
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n  n  n  n  n  n  n  n  n  d  n 


CHY  =  (2.0-AP*AP+2.0*AB-»AB) /AM/AM 

RT  =  SQRT(CHY»CHY-1 .0) 

CA  =  SQRT(ABS(CHY+RT) ) 

Ca  =  SQRT(ABS(CHY-RT) ) 

SS  =  A^fINl  (CA.CB) 

2AL  =  (2.0*ZETAA*ZETAB-»-(SS*SS*(ZETAA-2ETAB)-ZETAA-ZETAB) 

*  /CON  JG  ( ZETAA  )  )  /  (  5S*SS-<f  (  ZETAA-ZETAB  )  +ZETAA+ZETAB-2 . 0/ 

*  CONJG ( ZETAA ) ) 

ZBT  =  (  2 . 0»ZETAA^»-2ETAB-2AL*  (  ZETAA+ZETAB  )  )  /  (  ZETAA+ZE“ A3-2 . 0 

*  *ZAL) 

ZGM  =  SS*(ZETAA-ZBT) /(ZETAA-ZAL) 

415  (4RITE(G,245)  ZABST 

245  FORNAK //lOX, 'BEST  GUESS  FOR  ZETA  A  IS  ZABST  =  ',1P2E12.3) 
KRITE(S,24G)  ZBBST 

24G  FORMAT( //lOK,  'BEST  GUESS  FOR  ZETA  B  IS  ZBBST  =  ',lP2Ei2.3) 
KRITE(G,215)  ZETAA, ZETAB 

215  FORMAT(//  5X,  'ZETAA  =  ',1P2E13.5,'  ZETAB  =  ',2Ei3.5) 

WRITE (G,21B)  ZAL,ZBT,ZGM,SS 

21G  FORNATC//  5X,  'ALPHA  =  ',1P2E11.3,'  BETA  =  ',2E11.3, 

*  '  GAMMA  =  ',2E11.3,'  S  =  ',2E11.3) 

-INDING  THE  LOCATION  CF  BLADE-SURFACE  PC  I  NTS  IN  THE  ZETA  PLANE  C.'.L 
INUOLUES  PHI (THETA),  SINCE  R  =  1.  USE  SPLINE  INTERPOLATION 


PHKICS)  =  PHKl)  +  TPI 
I--  ( ICIRC.EQ.  1  )  GO  TO  211 
Ei 129)=e(l)+TPI 

CALL  CISPLN(PHI ,E,THT,F, 12S, 1 ,1,2) 
CALL  CISPLN(PHI ,E,X,F, 12S,2,KJMX,2) 
GO  TO  212 

211  CONTINUE 

£i 12S)=S(KJMX) 

CALL  CISPLN(PHI ,E,THT,F, 12S, 1 ,1,2) 
CALL  CISPLN (PHI,E,S,F,12S,2,KJMX,2) 

2 12  CONTINUE 
WRITE (S, 21 7) 


21.7  FORMATy/lOX,  'MAPPING  FROM  SMALL  OMEGA  -  PLANE  TO  ZETA  -  P,  OWF' 
♦  //  3X,  'KJ',1  IX,  'SMALL  OMEGA ',  23X ,' ZETA ', //T 


DO  54  K  =  1,KJMX 
R  =  EXP( Y(K ) ) 

ZX  =  CMPLX(R*COS(X(K ) ) ,R*SIN(X(K) ) ) 
ZVG  =  CMPLX(COS(F(K) ) ,SIN(F(K) ) ) 
IF(X.Ea. 1 .OR.K.EQ.KJMX)  ZYG  =  Z1 
ZnC(K)  =  ZYG 

54  WRITE(S,218)  K,ZX,ZYG 
218  FORMAT( 15, 1P4E15.5) 


zee  NOW  CONTAINS  ZETA  ON  THE  BLADE  SURFACES 
WRITE(G,202) 

NOW  DO  THE  MAPPING  FROM  THE  ETA  -  PLANE  TO  THE  KSI  TILDE  -  PLANE, 
WHERE  ETA/S  =  SN(KSI  TILDE) 


L-'L- ^  L -L 
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AK  =  SS-»SS 

AKQ  =  AX*AK 

AKP  =  SQRT( 1 .O-AKQ) 

AKM  =  AKP*AKP 

CALL  ELLPT(PB2,AK,RL,1) 

T3K  =  2.0*RL 
CTR  =  -RL 
CAPK  =  RL 

CALL  ELLPKPI  ,AKP,RL,  1  ) 

CAPKPM  =  RL 

kRITE(SF222)  AK,CAPX,AKPfCAPKPM 

222  F0RMAT(//10X,  'COMPLETE  ELLIPTIC  INTEGRALS  OF  K  AND  K  PRIMP' 

■S'  '  APf=’  pni  I  nLic«  /  .  XAJI-.  T 


^  '  ARE  AS  follows: ' , 

^  //lOX,  'K(  'fFIO.G,  '  )  =  ',F10.G,5X,  'K(  'fFIO.Bf  ')  =  'r-'OB'' 
DO  55  I  =  1  fKJMX  ■  ■  ■  ' 

ZAC  I)  =  ZGM-*(2CC<  I  )-2AL)  /  (ZCCC  I  )-ZBT) 


!A(I)  NOW  HOLDS  ETA 


=  ZACD/SS 
I  =  REALCZTD) 

■  =  AIMAGCZTD) 

I  =  TAU*TAU 
I  =  DLT»DLT 

■=1.0+  AKQ*(TSQ  +  DSQ) 

=  SQRTC ( 1 .0-AKQ*TSQ)*< 1 .0-AKQ*TSQ) 


AKQ»DSGi-»(2 .0» 


*  ( 1 .0+AKQ+TSQ) 

*  +AKQ*DSQ)) 

BRQ  =1.0+  TSG  +  DSQ 

BRT  =  SQRTC < 1 .0-TSQ)*( 1 .0-TSQ)  +  DSQ* C DSQ+2 . 0+2 . 0»TSQ ) ) 
A_M  =  CBRQ-BRT)+CART-RT) /4.0/AKQ/TSQ 
SGA  =  CTSQ  +  DSQ  -ALM) 

SGA  =  SGA/ CSGA+1 . O-ALM+AKQ+CTSQ+DSO) ) 

IF CTAU.EQ.0.0)  GO  TO  403 
RTALM  =  SQRTCALM)*TAU/A3SCTAU) 

GO  TO  404 

403  RTALM  =  0.0 

404  SGN  =  1,0 

IFCDLT.lt. 0.0)  SGIM  =  -1.0 
RTSGA  =  SGN^SQRTCSGA) 

RTALM  =  ASINCRTALM) 

RTSGA  =  ASINCRTSGA) 

CALL  ELLPTC-RTALMfAK  fRLfO) 

CALL  ELLPTCRTSGAfAKPfAGfO) 

I- CAG.GE.O.O)  go  to  57 
a:  =  -AG 
RL  =-RL  -  TBK 
57  ZAICI)  =  CMPLXCRLrAG) 

5:J  CONTINUE 


ZAICI)  NOW  HOLDS  KSI  HAT 


WRITECG.ZIS) 

‘19  F0RMATC////10Xf  'MAPPING  FROM  THE  ETA  -  PLANE  TO  THE  KSI  HAT  - 
+'  PLANE', 

*  ///  3Xr  'K J  '  , 15Xf  'ETA  '  f25X,  'KSI  HAT',//) 


»  V  >  ■ 


vv- 
*  *  ^ 
*  -•  V 
^ 


M 


"r* 

.•V'’  t 


.  -  ■ .  •  *  • 

•  */  ' 
'V-'. 
•vV.V.V 


DO  5G  K  =  1,KJMX 
kf?ITE(Sr218)  K  r2A(K  )  ,ZA1  (K  ) 

50  CONTINUE 
C 

C  NOW  SET  UP  A  GRID  IN  THE  KSI-HAT  PLANE,  AND  MAP  IT  BACK 
C  TO  THE  Z  -  plane: 

C 

ZA2(1)  =  ZTE 
2A2(KJMX)  =  ZTE 
ZA2(KJLE)  =  ZLE 
DO  91  KJ  =  2,KJLM 
K  *  KJLE  -  KJ 

51  ZA2(KJ)  =  ZP(K) 

DO  92  K  =  1 ,KJS 
KJ  =  KJLE  +  K 

92  ZA2(KJ)  =  ZS(K) 

C 

C  THE  ZA2  ARRAY  NOW  HOLDS  THE  BLADE-SURFACE  COORDINATES,  IN  THE  ORDER 
C  KJ  =  i:  TE 

C  KJ  =  2,KJLM:  pressure  side,  FROM  TE  TO  LE 

KJ  =  kjle:  le 

KJ  =  KJLP,KJMXM:  suction  side,  from  le  TO  TE 
KJ  =  KJMX:  TE  AGAIN 

KMXMl  =  KMX  -  1 
LMXMl  =  LMX  -  1 
HCPKPM  =  CAPKPM/2.0 
TWOCPK  =  2.0*CAPK 
THCPK  =  3.0*CAPK 
FCPK  *  4.0»CAPK 
ZAG  =  ZAL*2GM 
EKINU  =  1.0/EX 
SHXTE  =  REALCZAl ( 1 ) ) 

IF(SHXTE.GT.CTR)  SHXTE  =  SHXTE-FCPX 
5HK  =  -CAPKPM»CAPKPM/4.0/ (THCPK  +SHXTE) 

SHKINU  =  1.0/SHK 
IF( ISHEAR.EQ.O)  SHKINO  =  0.0 
DSHX  =  2.0/FL0AT(KMXMl ) 

DSHY  =  1 .0/FL0AT<LMXMl) 

ZIG  =  ZI*SG 
ZXYN  =  ZGAMMA*2N 
ZXYT  =  ZGAMMA*2T 
ZEN  =  CEXP(PI*ZXYN/SG) 

ZET  =  CEXP(Pr*ZXYT/SG) 

Zt- N1  =  l./ZEN 
ZETI  =  l./ZET 
SGP  =  SG/TPI 
SZ  =  0-5*SG 
ZEETE  =  ZA(1) 

ZLA  =  !.-»■  SS-*ZEETE 
ZLB  =  1 .  -  SS*ZEETE 
ZUC  =  ZEETE  +  SS 
ZLD  =  ZEETE  -  SS 

ZOQ  =  ZAM*(1./2LA  +  l./ZLB)#SS  +  ZAP* ( 1 . /ZLC  -  l./ZLD) 

ZOG  =  ZI* (SS*( 1 . /ZLA  -  1./ZL8)  -  l./ZLC  -  l./ZLD) 


w  '-J-^ '  I  ni  1  w 
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1 
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t; 


c 

c 

c 

c 


c 

c 


0  =  SG/TPI 
GG  =  -Q*ZDS/ZDG 

2P00T  =  4.»AK*(ZAM+AK»ZAP)-«-(ZAP+AK*ZAM) 

Zi^OOT  =  ZROOT-GG*GG*AKM*AKM/Q/Q 
RODT  =  REAL(ZROOT) 

ROOT  =  SQR'!' (  ROOT  ) 

ZSPB  =  (AKM*ZI*GG/Q  -  ROOT ) / ( 2 . *SS* ( ZAM+AK»ZAP ) > 

ALPO  =  ALP*180./PI 
GO  =  GG/Q 

WUP  =  SQRTC  1  .+GG*Ga-2.»GQ-»SIN(ALP>  ) 

KDN  =  SORTC 1 .+GG*GQ+2.*GQ*5IN(ALP) ) 

ULJP  =  CaS(ALP)/WUP 
UDN  =  UUP 

UDN  =  -(SIN(ALP)+GQ)/WUP 
UUP  =  -(SIN(ALP)-Ga)/WUP 
WR  =  WDN/WUP 

E2TA1  =  ATAN(UUP/UUP)*iaO./PI 
BETA2  =  ATAN(UDN/UDN)*iaO./PI 

'^R I TE  (  S  ,  22G  )  ALPO  ,  Q  ,  GG  r  ZSPB  ,  UUP  » UUP  ,  BETA  1  ,  UDN  ,  UDN  r  3ETA2  ,  WR 
22G  FORMAT( //lOK,  'CONSTANTS  FOR  I NCOMPRESS I BLE-FLON  SOLUTION  ARE' 
//:OXt  'ANGLE  OF  ATTACK  IN  THE  ETA  -  PLANE  ='rF10,5r 

*  //lOX.'Q  =  ',F10.5,'  GG  =  '.F10.5,'  ETA  <  STAG .  .  >  =  ', 

•ft  2F10.5,//5Xr  'INLET  U/WO . U/WO . 3ETA1  =',3F10.5, 

*  //5X,  'OUTLET  U/W0,U/W0rBETA2rW/W0  =  '  , 4F1 0 . 5 r / / ) 

A  =  1 

L  =  1 

WRITE{Gr202) 

kRITE(G,205) 

205  F0RMAT(/5X,  'MAPPING  OF  A  GRID  IN  THE  KSI-HAT  PLANE', 

ft  /SXr'AND  INCOMPRESSIBLE-FLOW  SOLUTION  IN  THE  X,V  PLANE', 
ft  //3X,'K  L',  8X,'KSI  HAT', 

ft  ISX,  'ETA ', IGX,  'ZETA ', 13X,  'SMALL  OMEGA', 

ft  10X,'Z  MAPPED ', 14X,  'ZXY  '  , 

ft  /13X,  'U/WO',GX,  'U/W0',7X,  ' PH I ' , 7X ,  '  PS I  '  , / / ) 

CALL  ZISPLN(ZSOMGA,ZCC,ZEE,ZFF,KJMX, 1 , 123, 1 ) 

K  -  LOOP  STARTS  HERE 
7G0  CONTINUE 

SHX  =  -1 .0+DSHX*FL0AT(K-l ) 

L  -  LOOP  STARTS  HERE 

770  CONTINUE 

SHY=DSHY*FL0AT<L-1 ) 

XIM=HCPKPM*( 1 .0-SHY) 

XIR  =  (XIM-HCPKPM)*ft2 

XIR  =  SHXTE+  TWOCPKft( 1 .OftSHX)ftXIRftSHKINU 
ZXI(L,K)  =  CMPLX(XIR,XIM) 
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SYPASS  IMAGE  CALCULATIONS  FOR  POINTS  THAT  FALL  ON  THE  IMAGES  Of 
PLUS  OR  MINUS  INFINITY  OR  THE  T.E. 


IF( ISHEAR.EO.O)  GO  TO  84 
IF (L.EQ. 1 .AND.K.EQ. 1 )  GO  TO  73 
IF(L.EQ. 1 .AND.K .EG.XMX)  GO  TC  73 
84  CONTINUE 

IFCL.EQ.LMX.AND.K.EQ. 1 )  GO  TO  74 
IFCL.EQ.LMX. AND.K .EQ.KMX)  GO  TO  74 
IF( lOE.EQ.O)  GO  TO  30 
IF(L. EQ.LMX.  AND.K  .EG, ;<MXH)  GO  TO  81 
GO  TO  80 


73  ZEETA(LrK)  =  ZEETE 

ZETA(L,X)  =  (ZBT»ZEETA(L,K )-ZAG)/(Z£ETA(L,K )-ZGM) 
ZOMA(L,K)  =  ZOMTE 
ZFNL(L.K)  =  CMPLXfSHX. SHY) 

ZXY(L-i<!  =  ZFE»ZGAMMA 
ZOEL<LrK)  =  ZBG 
ZFNP<L,K>  =  ZBG 
GO  TO  710 

74  ZEETA(LrK)  =  CMPLX ( SS r 0 . 0 ) 

ZETACLrK)  =  ZETAA 
ZOMA(LtK)  =  2MGA 
ZFNL<LrK)  -  CMPLX <SHX, SHY) 

ZOEL(LrK)  =  CMPLX (UDNrODN) 

ZFNP(L,K)  =  ZBG 

CO  TO  710 

81  ZEETA(L,K)  =  CMPLX < -SS , 0 . 0 ) 

ZFTA(LfK)  =  ZETAB 
ZOMA(L,K)  =  ZMGB 
ZFNL(L.K)  =  CMPLX (SHXr SHY) 

ZoEL(LrK)  =  CMPLX ( UUP, VUP) 

ZFNP<L,K)  =  ZBG 
GO  TO  710 


80  CONTINUE 

CALL  JCELFN(XIR,XIM, AKG, AKM,RLS,AGSr 1 ) 

ZEETA(L,K)  =  SS#CMPLX(RLS,AGS) 

IF(L.NE.LMX)  GO  TO  8G 

IF<K.LE.KMXH)  ZEETA(L,K)  =  ZEETA<LrK)  -  ZE?S 
I- (K. .GT.KMXH)  ZEETA(L,K)  =  ZEETA(L,K>  +  ZEPS 
83  CONTINUE 

Z£TA(L,K)  =  (ZBTttZEETA(L,K)-ZAG)/(ZEETA(L,K>-ZGM) 
IFfL.NE. 1 )  GO  TO  82 
ZEE( 1 )  =  ZETA(L,K ) 

CALL  ZISPLN(ZS0MGA,ZCC,ZEE,ZFF,KJMX,2, 1 , 1 ) 
ZGMA(L,K)  =  ZFF(l) 
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B2  IF  (ICIRC.Ea.O) 

CALL  GI*!ETA(A,B,ZOMA(L,K)  ,ZETA(L,K)  rZTAIMSRrES,  1  .OE-OOr  1  ) 

710  CONTINUE 
L  =  L  +  1 

IP(L„LE.LMX>  GO  TO  770 
K  =  K  +  1 
L  =  i' 

IF(K.LE.KN)<)  GO  TO  7S0 
IF( ICIRC.EO.O)  GO  TO  7B3 

C  AT  THIS  POINTr  GRID  CALCULATIONS  IN  THE  ZXI r  ZEETA,  AND  ZETA  PLAN 
C  ARE  COMPLETEr  FOR  ALL  L  AND  Kr  AND  IF  ICIRC  =  0,  THE  SMALL-GMEGA- 
C  PLANE  GRID  HAS  ALSO  BEEN  FOUND.  THE  FOLLOWING  CALL  CALCULA i ES  TH 
C  SMALL-ONEGA  GRID  FOR  THE  CASE  ICIRC  =  11 
C 

CALL  OMZTDRCDR.DI ,ZETABrZETA,ZOMArKMX,LMX,IOE) 

7B3  CONTINUE 
K  =  1 
L  =  1 

7G1  CONTINUE 
C 

C  the  :<  AMD  L  -  LOOPS  RESUME  HERE 

IF(L.EQ. 1 . AND.K .EG. 1)  GO  TO  711 
IF ( L. EG. 1 .AND. K .EG. KMX)  GO  TO  711 
I- (L.EG.LMX.AND.X .EG. 1 )  GO  TO  712 
IFIL.ED.LMX.AND.X.EG.KMX)  GO  TO  712 
Ii=<  I OE.  EG. 0)  GO  TO  7G2 
I^<L.EG.LMX.AND.K.EQ.KMXH)  GO  TO  712 
7E2  CONTINUE 

Z30M  =  2C-!KZ0MA(L,K  )-ZD)  /  (  ZOMA  ( L  » K  ) -ZB  ) 

RAD  =  CABS(ZBOM) 

ARG  =  ATAN2( AIMAu(ZBOM) rREAL(ZEOM) ) 

RADO  =  RAD**EXINV 
ARGO  =  EXINO*ARG 

ZBOMK  =  CMPLX(RADa#COS(ARGO) 7RADO*SIN(ARGO) ) 

C 

C  NOW  FIND  ZXY,  GIVEN  ZBOMK 

r 

ZXYlLrK)  =  SGP*CLOG(  (ZET-ZBOI'IK*ZEN)/(ZETI-ZBOMK*ZENI  )  : 

IFCL.NE. 1 )  GO  TO  33 

IF  (  (AIMAG(ZXY(L,K)  )-YY(K-l  ,L)  )  .LT.-SZ)  ZXY(L»K)  =  ZXY("_,k:  * 
IF( (AIMAG(ZXY(L,K ) )-YY(K-l ,L) ) .GT.SZ)  ZXY(L,K)  =  ZX'iL.K)  -ZI 
GO  TO  90 
33  CONTINUE 

IF(  (AIMAG(ZXY(L,K  )  )-YY(K  ,L-1  >  )  .LT.-SZ)  ZXY(LrK)  =  ZXYiL.K;'  +  Z 
IF(  (  AIMAG(ZXY(L,K  )  )-YY(K  ,L-1 )  )  .GT.SZ)  ZXY  (  L  ,  K  )  =  ZXY(L,;<)  -  2 1 
SO  CONTINUE 
72  CONTINUE 


ZFNL(L,K)  =  CMPLX(SHX.SHY) 
ZLA  =  1 .  +  SS#ZEETA(LrK ) 
ZLB  =  1.  -  SS*ZEETA(L,K ) 

7'  C  =  ZEETA  (  L,  K  )  +  SB 
ZLD  =  ZEETA(L,K)  -  SS 
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2ulE  =  Q*(SS*(  1  . /2LA+1  . /ZLB)*ZAM+ZAP-»(  1  ./ZLC-1  ./ZLD)  ) 
ZKE  =  ZWE+  ZI*GG*(SS*(1./ZLA-1./ZLB)-1./ZLC-1./ZLD) 
ZLA  =  CLOG(ZLA) 

Z^B  =  CLOG(ZLB) 

ZLC  =  C!_OG(ZLC) 

ZLD  =  CLOG(ZLD) 

X.C  =  REAL (ZLC) 

XLD  =  REAL (ZLD) 

YLC  =  AIMAG(ZLC) 

YLD  =  AIMAG(ZLD) 

IP(YLC.LT.O. )  YLC  =  YLC  +  TPI 
IFCYLD.LT.O. )  YLD  =  YLD  +  TPI 
ZLC  =  CMPLX(>(LCrYLC) 

D  =  CMPLX(XLD, YLD) 

ZFNP ( L , K )  =  Q* ( ZAM* ( ZLA-ZLB ) +ZAP* ( ZLC-ZLD ) ) 

ZfNP(L,K)  =  ZFNP(L,K )+ZI*GG*(ZLA+ZLB-ZLC-ZLD) 


C  CALCULATE  THE  UELOCITY 
C 

ZDR  =  ZDA 

ZDC  =  ZETA(L,K)  -  ZBT 
ZUC  =  ZGMXZAL-ZBT) /ZDC/ZDC 
7!!R  =  ZDR-^ZDC 
If-  (ICIRC.EQ.O) 

>  CALL  OMETA(A,B,ZOMA(L,K) ,ZETA(L,K) rZTANSRTBEr 1 . rZ) 

If-  ( ICIRC.EQ.  i  ) 

*CALL  QMDZETA(ZOMTErZl , DR . D I r ZOMA ( L , K ) rZETA(LrK ) -ZTAMSRr E5. 


2  OR 

s 

2DR/ZTANSR 

ZDE 

= 

20MA(L,K)  -  ZB 

2DE 

= 

ZC«’(ZD-ZB)  /ZDE/ZDE 

ZDR 

= 

ZDR/ZDE 

Z  = 

ZXY(LrX )/ZGA«MA 

Z2T 

= 

(Z-ZT) /ZDN 

Z7N 

= 

(Z-ZN)/ZDN 

ZT2 

= 

CSIN(ZPI»ZZT) 

Z  i  4 

= 

CSIN<ZPI*ZZN) 

2  :  1 

= 

CSIN(ZPI*(ZZN-ZZT) ) 

ZD- 

= 

ZPI*ZE0M*EX»ZT1/ZDN/ZT2/ 

ZDR 

= 

ZDR»ZDF 

ZUEL(L,K)  =  CONJG(ZwE*ZDR)/WUP 

GO 

TO 

711 

C 

C 

C 

712  Z)(Y(LrK)  =  2.0*ZA(Li'1X-l  )  -  ZA(LMX-2) 

7i  1  WRITE(3r224)  K  rLrZXI  (L,K  )  ,ZEETA(L.K)  ,2ETA(L,K  )  ,ZOMA(L,;<  )  r 
*  ZFML(LrK ) ,2XY(LrK ) 

U;"'  TE(S.225)Zi.''EL(L,K)  ,ZFi'IP(LrK) 

225  ~0RMAT(8X, 12F10.5) 

2P:L)  =  ZXY(LrK) 

X:<(K,L)  =  REAL(ZXY(LrK  )  ) 

YV(K,L)  =  AIHAG(ZXY(L,K )  ) 

7;  CONTINUE 
L  =  L  +  1 

I(- (L.LE.LNX)  GO  TO  7G1 
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|.:'?ITE(G,202) 

IF(PNCHZA)  WRITE(7,210)  ( 2A  ( L  >  .- L=  1  r  LMX  > 

K  =  K  +  1 
L  =  1 

IF(K.LE.KMX)  GQ  TO  7G1 
70  CONTINUE 

STOP 

100  F0RMAT(BF10.4> 

202  FORMAT (///) 

210  FORMAT( 1P4E20. 13) 

2:';4  FaRMAT(2I4, 12Fi0.5) 

END 


SUBROUTINE  C I SPUN ( Y , X , E , F , NP , IRTN , NRTN , NPD ) 

C 

C  THIS  SUBROUTINE  FITS  A  CUBIC  SPLINE  TO  A  FUNCTION  Y(X),  DEFINED  BY 
C  NP  PAIRS  OF  POINTS.  THE  SECOND  DERIVATIVE  OF  THE  FUNCTTCN 
C  IS  !=ERIODIC.  IF  NPD  =  1,  THE  FUNCTION  ITSELF  IS  ALSO 
C  PERIODIC.  WHILE  IF  NPD  =  2,  THE  FUNCTION  INCREASES  BY  2.-»-PI 
C  EVERY  PERIOD,  IN  WHICH  CASE  THE  DATA  PASSED  TO  THIS 
C  SUBROUTINE  SHOULD  HAVE  THE  PROPERTY  THAT  Y(NP)  =  Y(l)  *  2.*Pi; 

C  THE  SUBROUTINE  WILL  THEN  SET  : 

C  Y(NP+1)  =  Y(2)  +  (NPD-i  )*2.-*PI ,  AND 

C  X(MP+i)  =  X(NP)  +  X)2)  -  X(l) 

C  SOLUTIONS  FOLLOW  PAGES  9  -  15  OF 

C  THE  THEORY  OF  SPLINES  AND  THEIR  APPLICATIONS,  BY  J.  H.  AHL3ERG, 

C  E.  N.  NILSDN,  AND  J.  L.  WALSH,  ACADEMIC  PRESS,  1SB7 

C 

C 

C  NOTE  THAT  Y,X,E,AND  F  ARE,  RESPECTIVELY,  ORDINATE,  ABSC I S3A , ABSC I SSA , 
C  AND  ORDINATE. 

C 

IMPLICIT  REAL(A-H,0-Y) ,CCMPLEX(Z) 

DIMENSION  Y( IGO) ,X( IGO) ,E< ISO) ,F( ISO) ,BDA( IBO) ,EM(1G0) ,H( IGO) 

D ! MENS ION  S  (  1  GO ) , T ( 1 GO ) , V ( 1  GO ) , D ( 1  GO ) 

D^TA  TPI/G.2831S53071795/ 

C 

C  THIS  SECTION  (ENTERED  WHEN  IRTN  =  1)  USES  THE  NP  PAIRS  OF  INPUT 
C  COORDINATES  X  AND  Y  TO  FIND  THE  COEFFICIENTS  OF  THE  SPLINE  FIT. 

C 

C  THESE  COEFFICIENTS  -  HERE  CALLED  EM(KJ)  -  ARE  THE  SECOND  DERI  VAT  I '-'Eb 
C  OF  THE  FUNCTION. 

C 

f ; 

IF( IRTN. EG. 2)  GO  TO  20 
NPM  =  NP  -  1 
N  =  NP  +  1 
DO  1  KJ  =  2,NP 

1  H(KJ)  =  X(K J)  -  X(KJ-1 ) 

H  (  N  )  =  H  (  2  ) 

DO  2  KJ  =  2,NP 

2  BDA(KJ)  =  H(K J+1 > / (H<K J ) +H(K JTI ) ) 
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Ed)  =  0.0 
F  ;  1 ;  =  0 . 0 
S  (  1  )  =  i  .  0 
Dl'  3  KJ  =  2»NP!»i 

DM  =  2.0  +  a.O  -  BDA(KJ)  )»E(KJ-1  ) 

E(KJ)  =  -3DA(KJ)/DN 
D  ■  '  '<  J  )  = 

*  S.0*( (Y(KJ+1 )-Y(KJ) )/H(KJ+l ) 

»  -( Y(KJ )-Y(K J-1) ) /H(K J ) ) /(H(K J)+H(K J+i ) ) 

F.!  KJ)  =  -(  1 .0-8DA(KJ)  )*S(KJ~1)/DIM 

3  -(KJ)  =  (D(K  J  )-(  1 .0-3DA(K  J  )  )»F(KJ-1  )  ) /DIM 
Y  (  M  )  =  Y  (  2 ) 

IF(NPD.EQ.2)  Y(N)  =  Y(N)  +  TPI 
D(MP)  =  B.0*( (Y(N)-Y(NP) )/H(2) 

»  -(Y(NP)-Y(NPM) )/H(NP) )/(H<NP)+H(2) ) 

N^MM  =  NP  -  2 
T(NP)  =  1.0 
'J(.MP)  =  0.0 
DO  B  I  =  l.NPMM 
KJ  =  NP  -  I 

T(KJ)  =  E(KJ )*T(K J+1 )  +  S<KJ) 

U(KJ)  =  E(KJ)*0(KJ+1 )  +  F(KJ) 

G  CONTINUE 

F^(NP)  =  (D(NP)-BDA(NP)»0<2)-( 1 .O-BDA(NP) )»0(NPr ) ) / 
fr  (BDA(NP)*T(2)+( 1 .O-BDA(NP) )#T(N?M)  +  2.0) 

DO  4  I  =  i,NPM 
KJ  -  N?  -  I 

4  EM(KJ)  =  E(KJ)#EM(KJ+1 )  +  F(KJ)  +  S(KJ)»EM(NP) 

5  .ZE  =  ABS()<(NP)-X(1  )  )/10.0 

rytTL'RN 

C 

C  THIS  SECTION  (ENTERED  WHEN  IRTN  =  2)  RETURNS  NRTN  INTERPOLATED 
C  0ALUE5  OF  THE  ORDINATE  F  AT  ASSIGNED  VALUES  OF  THE  ABSCISSA  E. 

C 

20  KJ  =  2 

DO  21  J  =  1,NRTN 
=  E<J) 

24  IF(A.LE.X(KJ) )  GO  TO  23 
K  =  K  J  +  1 

IFCKJ.LE.NP  )  GO  TO  24 
DE  =  A  -  )<(NP) 

IF(DF.GT.SIZE)  WRITE(Gr200)  J , E ( J ) , NP , X ( NP ) 

200  FORMAT < //lOX,  'WARNING  -  ENTRY  IN  CI5PLN  EXCEEDS  END  OF  BASE 

^  'ARRAY /5X,  'E(  13,  '  )  =  'dPElB.Sr'  EXCEEDS  X(',I3,')  =  ' 

*  EIB.B) 

D-  =  A  -  X (  1  ) 

IF(DF.LT. (-SIZE) )  WR I TE ( B , 20 1 ) J , E ( J ) , X ( 1 > 

201  FORMAK // lOX, 'WARNING  -  ENTRY  IN  CISPLN  IS  LESS  THAN  ^HE  FIRST', 

*  '  BASE  POINT', 

*  /5X, 'E( ' , 13, ' )  =  ',iPElG.8,'  13  LESS  THAN  X(l)  =  ', EIB.B) 

KJ  =  NP 

23  DXA  =  X(KJ)  -  A 
DXB  =  A  -  X(KJ-l) 

CBA  =  DXA*DXA*DXA 
CBB  =  DXB*DXB*DXB 

F(J)  =  (EM(KJ-1)*CBA+EM(KJ)*CBB)/B.0/H(KJ) 

*  +  DXA*(Y(KJ-1)-EM(KJ-1)*H(KJ)*H(KJ)/S.0)/H(KJ) 

*  +  DXB#(Y(KJ)-EM(KJ)*H(KJ)*H(KJ)/B.O)/H(KJ) 

21  CONTINUE 
RETURN 
END 


V  V  V 

V.-.'C'-Lv 

VV,\  -v- 


»  M  •  - 
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SUBROUTINE  ELLPT ( AR , AK . ANS r KOMP ) 

IMPLICIT  REAl(A-H,0-Y) ,C0MPLEX(2) 

D;  PENSION  SQ(12) 

DATA  K./l/ 

DATA  PI/3. 1415926535897/ 

C 

C  THIS  SUBROUTINE  EVALUATES  THE  ELLIPTIC  INTEGRAL  OF  THE  "IRST  KT'-iD, 
WITH  ARGUMENT  AR  (AN  ANGLE  IN  RADIANS),  AND  PARAMETER  AK  (A  REAL 
C  NUMBER). 

C  THIS  EVALUATION  USES  EG.  (14)  OF!  Y.  L.  LUKE,  -'APPROXIMATIONS 
C  FOR  ELLIPTIC  INTEGRALS',  MATH.  COMP.,  MOL.  22  (JULY  ISSS),  PP  327- 
C  G34,  WITH  N  =  12. 

C  KOMP  =0,1  FOR  THE  INCOMPLETE,  COMPLETE  INTEGRAL,  RESPEC'IMELV , 

C 

IF(K.GT.l)  GO  TO  11 


TNP  =  25.0 

DO  10  M  =  1,12 

THM  =  PI#FLOAT(M)/TNP 

S  =  SIN(THM) 

10  5G(M)  =  S-»-S 
1  i  AKK  =  AK»AK 
S'-  =  0.0 

IF(KCM?.£Q. i )  GO  TO  40 
Tn  =  TAi\i(AR) 

D-.  20  M  =  i,12 

SO  =  SQRT(  1 .0-AKK«-SQ(M)  ) 

T  =  ATAN(SG*TN) 

20  SM  =  SM  -  T/SG 

AN3  =  (AR  4-  2.0*SM) /TNP 
RETURN 

40  DO  41  M  =  1,12 

SG  =  SQRT( 1 .0-AKK*SQ(M) ) 

41  SM  =  SM  +  1,0/SG 

ANS  =  PI*( 1 .0+2.0*SM) /2.0/TNP 

RETURN 

END 


^  SUBROUTINE  JCELFN ( RL , AG , AKQ , AKM ,RLS , AGS , M ) 

'  •  WHEN  M . EG , 1 , 

i.  THIS  SUBROUTINE  RETURNS  THE  JACOBIAN  ELLIPTIC  SINE 
:•  OF  A  COMPLEX  ARGUMENT: 

!^1-5  +  I*AGS  =  SN(RL  +  I*AG,K) 

C  USING  THE  ARITHMETIC  -  GEOMETRIC  MEAN  FORMULA  (SEE  P  '^71,  HAMPRnr.K 
C  Oh  MA.HEMATICAL  FUNCTIONS,  ED.  BY  M.  ABRAMOWITZ  AND  A  StSuN 
r  BURtAU  OF  STANDARDS,  APPLIED  MATHEMATICS  5£RI-5^  5S, 

r  ADDITION  FORMULA  FOR  THE  SN  (SEE  FQUATtgn  •’-'t^O' 

ELLIPTIC  INTEGRALS  FOR  ENGINEERS  AND  ‘ 

i-  AKP  -  SPRINGER  MERLAG,  1S54). 

^  quantities  returned  ARE  THE  REAL  AND  IMAGINARY  PARTS 

U  OF  THE  PRODUCT  CN  (  .  .  )  #DN  (  .  .  )  ,  WHICH  IS  THE  DERIMATIME  DF  THE  5N>'.,) 


■'i'.'R.-: 


•TlV.vS' 
►V-  J- 

I-'  •  *  * 


*  'N.*  •. 

o.  •**, 


IMPLICIT  REAL(A-H,0-Y) .CQMPLEXCZ) 

DIMENSION  A(20) ,3(20) ,C(20) »PH(20> 

K  =  1 

A( 1 )  =  1.0 
B(l)  =  SQRT(AKM) 

C(l)  =  SQRT(AKQ) 

5  00  e  I  =  2,20 

A(I)  =  ( A( I-l )+B( I-l ) ) /2.0 
B!I)  =SQRT(A( I-1)*B( I-l ) ) 

CiI)  =  (A( I-l)-B( I-l ) ) /2.0 
IP (ABS(C( I ) ) .LT. 1 .OE-OS)  GO  TO  7 
S  CONTINUE 

WRITE (6,200)  RL,AG 

200  FORMAT( ///lOX, ' JCELFN  FAILED  TO  CONVERGE  FOR  Z  =  ',lP2Ei5,d) 
S  I  OP 

7  NM  =  I-l 
N  =  I 

IF (K. EG. 2)  GO  TO  20 
Ph(N)  =  A(N)*RL*2**NM 

15  DO  11  L  =  1,NM 
J  =  N  -  L 

II  PH(J)  =  (  PH(  J+1  )+ASIN(C(  J-t-1  )*SIN(  PH(  J  +  1  )  ) /A(  J+1  )  )  ) /2.0 
IP- (K.  EG.  2)  GO  TO  40 
SA  <  =  SIN ( PH ( 1  )  ) 

CNK  =  COS( PH( 1 ) ) 

DNK  =  CNK/COS( PH(2)-PH( 1 )  ) 

A  =  2 

■^MP  =  B(l) 

8(1)  =  C ( 1  ) 

C(  1 )  =  TMP 
GO  TO  5 

20  PH(N)  =  A(N)*AG*2**NM 
Gil  TO  15 

40  SNP  =  SIN( PH(  1  )  ) 

CNP  =  C0S(PH(1 ) ) 

DNP  =  CNP/C0S(PH(2)-PH( 1 ) ) 

D\M  =  1.0  -  SNP*SNP*DNK*DNK 

IF (M. EG. 2)  GO  TO  50 

RlS  =  SNK*DNP/DNM 

AGS  =  CNK*DNK*SNP*CNP/DNM 

RETURN 

50  RLA  =  CNK*CNP 

AGA  =  SNK*DNK*SNP*DNP 

Ri  B  =  DiMK«-CNP*DNP 

ATS  =  AKO*SNK*CNK*SNP 

RLS  =  (RLA*RLB-AGA*AG8) /DNM/DNM 

AGS  =  (-RLA*AG3-RLB*AGA)/DNM/DNM 

Rl  ’‘URN 

ErcO 
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SUBROUTINE  OMETA ( A , B , ZMGA  f  ZTAGS . ZTANSR  r  N . EPS  f  M ) 

}  ~  -'LICIT  REAL(A-H,0-Y)  ,COMPLEX(Z) 

Di' MENS  I  ON  A(G5)  FB<e5)  fZC(B5) 

IF  M.EQ.Of 

THIS  SUBROUTINE  USES  NEwTON  -  RAPHSON  TO  FIND  ZETA(OMEGA):  Zf-GA  IS  A 
KNOWN  VALUE  OF  OMEGA f  ZTAGS  IS  THE  INITIAL  GUESS  AT  ZETA f  ZTANSR  IS 
THE  SOLUTIONf  AND  EPS  IS  THE  TOLERANCE  ON  THE  ANSWER. 

M  IS  RESET  TO  5  AS  AN  ERROR  RETURN  IF  THESE  ITERATIONS 
DO  NOT  CONVERGE. 

IF  M.EQ.I.f  this  SUBROUTINE  RETURNS  THE  VALUE  OF  OMEGA  (IN  ZMGA) 

FOR  A  GIVEN  VALUE  OF  ZETA  (IN  ZTAGS). 

IF  M.EQ.Zf  the  quantity  returned  (IN  ZTANSR)  IS  D  OMEGA/D  ZETAf  FOR 
GIVEN  VALUES  OF  OMEGA  (IN  ZMGA)  AND  ZETA  (IN  ZTAGS) 

Z1  =  CMPLX( 1 .OfO.O) 

ZETA  =  ZTAGS 
NM  =  N  -  1 
IT  =  1 
5  CONTINUE 

ZSMA  =  CMPLX ( A ( N ) , B ( N ) ) 

I~(M.EQ. 1 )  GO  TO  22 
ZSM8  =  ZSMA-»FLOAT(NM) 

If(M.EQ.2)  GO  TO  40 
DO  10  J  =  1,NM 

ZSMA  =  ZETA»ZSMA  -f  CMPLX  ( A  (  N- J  )  ,  B  ( N- J  )  ) 

ZsHB  =  ZETA+ZSMB  +  CMPLX ( A ( N-J ) f B ( N-J )) ♦FLOAT ( MM- J ) 

10  CONTINUE 

ZEXP  =  CEXP(ZSMA) 

ZG  =  ZETA^ZEXP  -  ZMGA 
ZOG  =  ZEXP*(Z1  +  ZSM8) 

ZFTOLD  =  ZETA 
ZETA  =  ZETOLD  -  ZG/ZDG 
IPCCABSCZETA-ZETOLD) .LT.EPS)  GO  TO  20 
IT  =  IT  +  1 

I~(CABS(ZETA) .GT. 1 .0)  ZETA  =  0.9  *ZETA/CABS ( ZETA > 

IF( IT.LE.50)  GO  TO  5 
M  =  5 
RETURN 

20  ZTANSR  =  ZETA 
RETURN 

22  DO  21  J  =  IfNM 

ZSMA  =  ZETA»ZSMA  +  CMPLX ( A ( N- J ) f B ( N- J ) ) 

CONTINUE 

ZtXP  =  CEXP(Z5MA) 

ZMGA  =  ZETA*ZEXP 
RETURN 

40  DO  41  J  =  IfNM 

41  ZSMB  =  ZETA*ZSMB  +  CMPLX ( A ( N- J ) f B ( N- J > ) *FLOAT ( NM- J ) 


CMPLX(A(N-J) fB(N-J) ) 


CMPLX(A(N-J) fB(N-J> )*FL0AT(NM-J ) 


Z ! ANSR 
RETURN 
END 


ZMGA* ( Z 1 +ZSMB ) /ZTAGS 


.-v 


r,  •  -  «  ,  -r  ,  ^  ' 
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SUBROUT I NE  SHAPE (CfH.G.EX,ZP»ZS,KJSrKJP) 
IK^LICIT  REAL(A-H,0-Y) ,COMPLEX(Z) 
DIMENSION  ZS(80) rZP(80) 

DO  10  K  =  IxKJS 
10  READ (5, 100)  ZS(K) 

Df  20  K  =  1,KJP 
20  READ(5rl00)  2P(K) 
iriO  FORMAT (8F1 0.0) 

RETURN 

END 


SUBROUTINE  SHUFL ( N » ZA , ZCC ) 

C 

C  THIS  SUBROUTINE  TAKES  THE  N  COMPLEX  UALUES  IN  ARRAY  2A ,  WHICH  WERE 
C  COMPUTED  AND  STORED  IN  REVERSE  BINARY  ORDER  BY  "FT  AND  ''SHUF!=-L 
C  THEM  INTO  PROPER  ORDER  USING  ARRAY  ZCC  FOR  INTERMEDIATE  STORAG 
C  N  IS  ASSUMED  TO  HAVE  THE  FORM  2**M. 

C 

IMPLICIT  REAL(A-H,0-Y) ,C0MPLEX(2) 

Dimension  za(G5) ,zcc<b5) r ial<6) fKR(S4) 

D'^'TA  KALL/0/ 

DATA  IAL/G*0/ 

IF  (KALL.EQ.l)  GO  TO  10 
KA_L  =1 
D'  341  JP  =  1,N 
J  =  JP  -  1 
lAL(G)  =  J/32 
J  =  J  -  32*IAL(G) 

IAL<5)  =  J/IG 
J  =  J  -  1G*IAL(5) 

I^.v_(4)  =  J/S 
J  =  J  -  8*IAL(4) 

IAL(3)  =  J/4 
J  =  J  -  4*IAL(3) 

IAL(2)  =  J/2 
J  =  J  -  2*IAL(2) 
lAL(i)  =  J 

341  KRUP)  = 

32*IAL<  1  )  +  1G*IAL(2)+8#IAL(3)+4*IAL(4)+2*IAL(5)+IAL(G:' 

10  DO  342  J  =  1,N 

342  2CC( J)  =  ZA(KR( J)+l > 

DO  3B0  JP  =  IrN 

3I-.0  2A(JP)  =  ZCC(JP) 

RETURN 
END 
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SUBROUTINE  THDGRK 


C  THIS  SUBROUTINE  MAPS  AN  OVAL  TO  A  UNIT  CIRCLE,  USING  A  VARIANT  OF 
C  THE  THEODORSEN-GARRICK  TRANSFORMATION  AND  FAST  FOURIER  TRANSFORM 
C  TECHNIQUES.  (SEE  REFERENCES  AT  BEGINNING  OF  MAIN  PROGRAM), 
r 

IMPLICIT  REAL(A-H,0-Y) ,COMPLEX(Z) 

COMMON/TGINTG/N.NPl ,NP2.N2,NB2,NB2P1 , IP, IPMX, ITP , I WK , I MX , K JMX 
COMMQN/TGCMPX/Zl ,ZW2N,ZI ,ZNN,ZA,ZCC,ZA1 ,ZA2 

COMMON/TGDBLE/PBN,OM,OMM,ANGERR,A,B,E,F,THT, PHI ,X, Y,BETA,S,DR, DI 
DIMENSION  X< ISO) , Y( ISO) ,E( ISO) ,F( ISO) ,THT( ISO) , 

*  3ETA( ISO) ,S( ISO) 

D I MENS ION  PH I ( 1 30 ) , A ( S5 ) , B ( G5 ) , ZCC ( G5 ) , DR ( B5 ) , D I ( B5  )  , 

*  ZAC  IBS) ,2Ai (IBS) ,ZA2(1S5) , 

*  ■  IWK(7) 

DIMENSION  ITP(IOO) 

C 

301  DO  300  I  =  1,N2 
300  PHI ( I ) =PBN*FLOAT( I-l ) 

C 

DO  10  I=1,S5 
ZCCCI)  =  0. 

A( I )  =  0. 

B  (  I  )  =  0 . 

10  CONTINUE 

DO  11  1=1, ISO 

11  Ed)  =  0. 


WRITE(G,401 ) 


1433? 


WRITE(G,401 ) 

401  FORMAT! //lOX, 'PROGRESS  OF  PHI  /  THETA  ITERATIONS  IS  AS  FOLLONS:'// 
1  3X,  'IT' ,4X,  'DEMX' ,4X,  'NO.  OF  THETA  REVERSALS') 

C 

IT  =  1 
C 

C  FIRST  GUESS  IS  THETA  -  THETA ( TRAILING  EDGE)  =  PHI 

C 

DO  3 1 S  K  =  1  ,  N2 
3  15  E(K)  =  X(l)  +  PHKK) 

305  DEMX  =  0.0 
C 

C  USE  AJ  AND  BJ  TO  GET  NEXT  APPROXIMATION  TO  THETA 


ZCCCI )  =  CMPLX(0.0,0.0) 

DO  302  JP  =  2,N 

302  ZCC(JP)  =  CMPLX(B( JP)/2.0,-A( JP)/2.0) 
ZCC(NPi)  =  CMPLX(0.0,0.0) 

ZW  =  Z1/ZW2N 
DG  303  NP  =  1,NB2P1 

ZAKNP)  =  ZCC(NP)  +  CDNJG(ZCC(NP2-NP)  ) 
Zul  =  2W*2W2N 
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C 

c 

c 

c 

c 

c 


303  ZA2(NP)  =  ZW*(ZCC(NP,'-C0NJG(ZCC(NP2-NP)  )  > 

DO  304  IMP  =  1,NB2P1 

304  ZA(NP)  =  ZAKNP)  +  ZI*ZA2(NP) 

DO  309  NP  =  2,NB2 

305  ZA(NP2-NP>  = 

*  C0NJG(ZA1 (NP) )  +  ZI*C0NJG(ZA2(NP) > 

CALL  FFKZA.S,  IWK)- 

CAlL  SHUFL(N,ZA,ZCC) 

BU)  =  X(l)  -  PHI(l)  -  REAL(ZA(i)) 

DiJ  30G  K  =  lrG4 
TIMP  =  E(2*K-1) 

E:2*K-1)  =  REAL(ZA(K))  +  PHI<2*K-1)  +  B(l) 

ThT(2*K-1)  =  E(2*K-1) 

DEN  =  AB3(TMP-E(2»K-1 ) ) 

IF(DEM.GT.DEMX)  DEMX  =  DEM 
E(2*K-1)  =  □M*E(2»K-1)  +  0I1M*TMP 
Ti’iP  =  E(2*K) 

E(2*K)  =  AIMAG(ZA(K))  +  PHI(2*K)  +  B(l) 

THT(2*K)  =  E(2»K) 

D'N  =  ABS<TMP-E(2*K ) ) 

I- (DEN. GT. DEMX)  DEMX  =  DEM 
EiZ^K)  =  0M»E(2*K)  +  QMM-frTMP 
30G  CONTINUE 

IF( IT.NE. ITP( IP) )  GO  TO  940 
IP  =  IP  +  1 
iA;--'ITE(G,a41  )  IT 
841  FORMAT (/3Xr 

*  'THETA  BEFORE  AND  AFTER  RELAXATION  AT  IT  =  ',14) 
W?ITE(G,B32) (THT( I ) , 1=1 , 128) 

*  WRITE(G,a32) (E( I ) , 1=1 , 128) 

S32  FORMAT( 1P10E13.5) 

8-0  CONTINUE 

IE ( IT.GE. ITP( IPMX) )  STOP 

NOW  USE  THESE  THETAS  TO  GET  THE  NEXT  APPROXIMATION  TO  lN  R i K ) 

CALL  CISPLN(Y,X,E,F,KJMX,2, 128, 1 ) 

NOW  FIND  THE  AJ  AND  BJ  COEFFICIENTS  CORRESPONDING  TO  THE  LN  R(K.  DATA 
DO  310  JP  =  1 ,N 

310  ZA(JP)  =  CMPLX(F(2*JP-1 ) ,F(2*JP) ) 

DO  307  NP  =  1,N 

307  ZA(NP)  =  CONJG(ZA(NP) ) 

CALL  FFT(ZArG, IWK ) 

DO  308  JP  =  1,N 

3'>8  ZA(JP)  =  CONJG(ZA(  JP)  )/ZNN 
CA_L  SHUFL(N,ZA,ZCC) 

ZA(G5)  =  ZA(1) 

DO  31 1  NP  =  1 ,NB2P1 

ZAKNP)  =  (C0NJG(ZA(NP2-NP)  )  +  ZA(NP))/2.0 

311  ZA2(NP)  =  ZI*(C0NJG(ZA(NP2-NP) )  -  ZA(NP))/2.0 
ZU!  =  2W2N 

DO  312  NP  =  1 ,NB2P1 
ZW  =  ZW/ZW2N 
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3x2  ZCC(NP)  =  (ZAi(iviP)fZA2(NP)*ZW)/2.0 
Zv  =  Zl/ZW 
DC)  313  I  =  2,NB2 
2ui  =  ZW*ZW2N 
=  N32  +  I 

ZS3  =  (ZAl  (NP2-.NP)+ZA2iNP2-NP)*2W)/2.0 
3' 3  ZCC(NiP)  =  C0NJG(ZB8) 

A'l)  =  REALCZCCd)) 

ZCC(l\iPl  )  =0.5  *C0NJG(ZA1  <  1)-ZA2(  1  >  ) 

A(S5>  =  REAL ( zee (S5) ) 

DP  314  NP  =  2,N 
A ( NP )  =  2 . 0*REAL ( Zee ( NP ) ) 

31.4  B(NP)  =  2.0*AIMAG(ZCC(NP)  ) 

CHEeK  FOR  CONOERGENeE 

IF(IT.EQ.l)  DEMX  =  1.0 
IMRO  =  0 

DO  809  LL  =  2,128 

B09  IF(E(LL) .LE.E(LL-1 ) )  NRU  =  NRO  +  1 
kP ! TE ( e , 400 )  IT, DEMX , NRO 
400  FQRMAT( 15, 1PE12.3, 18) 

IF(DEMX.LT. ANGERR  )  GO  TO  31 B 
IT  =  IT  +  1 

IF( IT.LE. IMX)  GO  TO  305 
k"'ITE(6,213) 

213  F0RMAT(//5X,  'ITERATIONS  FOR  A  AND  B  DID  NOT  ODNOERGE '  ) 

S  :'0P 

315  WRITE (5,203)  IT,OM,DEMX 

203  F0RMAT(//10X, 'A  AND  B  ITERATIONS  CONVERGED  AT  IT  =  ',13, 
«  '  USING  OM  =  ',FG.3,'  .  MAXIMUM  ANGULAR  ERROR  =', 

iPE12.3,'  RADIANS') 

WRITE(G,220) 

220  FORMAT(//  4X  ,' l< 1 1 X  ,' THETA ',  15X ,' PHI ',// ) 

I/O  SS  K  =  1,128 

SO  WRITE(G,22i)  K , E ( K ) , PH I ( K  ) 

221  FORMATdS,  1P2E20.G) 

RETURN 

END 


SUBROUTINE  BAUGRK < RMAX , RM I N ) 

THIS  SUBROUTINE  MAPS  AN  OVAL  TO  A  UNIT  CIRCLE,  USING  THE 
HALSEY- JAMES  AND  FAST  FOURIER  TRANSFORM  TECHNIQUES. 

IMPLICIT  REAL< A-H,a-Y) ,COMPLEX(Z) 

CnMMON/TGINTG/N,NPl , NP2 , N2 , NBZ , NB2P1 , IP, IPMX, ITP, INK , IMX,K JMX 
COMMON/TGCMPX/Zl ,ZW2N,ZI , ZNN , ZA , ZCC , ZA 1 , ZA2 

COMMON /TGDBLE/PBN , OM , OMM , ANGERR , A , B , E , F , THT , PH  I , X , Y , BE  FA , S , DR , 

Dimension  x( iso) ,y(iso) ,E( ibo) ,F(ieo) ,tht( iso) , 

*  BETA ( ISO) ,S( ISO) 


DIMENSION  PHI ( 130) ,A<G5) ,B(G5) ,2CC(G5) ,DR(S5) rDI (G5)  r 
*  ZA(1G5) rZAi ( 1B5) ,ZA2( 1G5) , IWK (7) 

d:mension  IIP (100) 

DO  8  I=lrN2+l 
8  PHI ( I ) =PBN*FLOAT( I-l  ) 

DO  9  I=lrNPl 

zc,c(  I )  =cmpl;<(0.  ,0. ) 

D ' ( I ) =0. 

D:''(  I  )=0. 

S  CONTINUE 

DO  10  I =1,1 GO 

10  E( I )=0. 

IT  =  0 

FIRST  GUESS  IS  ABS ( D0KEGA/D2ETA )  =  ( RMIN+RMAX ) /2 . 

RAO=(RMAX+RMIN) /2. 

DO  11  J=1,N2+1 

11  F(J)=RAM 
HPBN=PBN/2. 

WRITE(G,401 ) 

401  FORMATC //lOX,  'P'ROGRESS  OF  PHI/S  ITERATIONS  IS  AS  FOLLOWS 
^3X,  '  IT  '  ,4)<,  'DSMX  '  ,GX,  'NO.  OF  ARG  ( DOMEGA/DZETA  )  REOERSALS 
S99  CONTINUE 
:-;=iT-t-i 

IF  (IT.LE. IMX)  GO  TO  801 
;*.:PITF(G,800) 

S.-.'O  FORMATi //5X,  '  ITERATIONS  FOR  DI  AND  DR  DID  NOT  CONOERGE '  ) 
S'!  OP 

SOI  CONTINUE 
DSMX=0. 

USE  ABS (DOMEGA/DZETA)  TO  APPROXIMATE  THE  NEXT  5 
E(  1  )=0. 

DO  12  I=2,N2+1 
THT(I)=E(I) 

E( I ) =E( I-l )+HPBN*(F( I )+F(  I-l )  ) 

DL  13  I=2,N2+1 

E(  I  )  =  S(KJMX)*E(  I  )  /  E(N2-^1  ) 

DS  =  A3S ( tut ( I ) -E (  I  )  ) 

I-  (DS.LE.DSMX)  GO  TO  15 
ESMAX=E( I ) 

D'^MX  =  DS 
I  ;'KMAX=I 
15  CONTINUE 
-'■',P  =  THT  (  I  ) 

THT( I ) =E( I  ) 

Ell) =0M»£ ( I ) +OMM*TMP 
13  CONTINUE 

IF  ( IT.NE. ITP( IP) )  GO  TO  840 
TFT! 1 ) =E( 1 ) 


t-jRITE  (8-841  )  IT 

841  FORMAK 13X- 'S  BEFORE  AND  AFTER  RELAXATION  AT  IT  = 
i'(RITE(B-332)  ("HT(  I )  - 1  =  1  ,i\|2+l ) 

K';lTE(S-832)  (E(  I)  ,  1  =  1  ,N2+1  ) 

632  FORNAT( 1  PI OE 13. 5) 

640  CONTINUE 

IF  ( IT.GE. ITP( IPMX) )  STOP 

CHECK  FOR  CONVERGENCE 

IF  (IT.EQ.i)  DSMX=1. 

NRV  =  0 

DO  809  LL=2-N2 

809  IF  (E(LL)  .LE.E<LL-1  )  )  NRV  =  NRV^-1 
WR : TE ( G -  400 )  IT, DSMX , NRV 
4  00  FORNATdS,  1PE12.3,  IB) 

IF  (DSNX.LT.ANGERR)  GO  TO  SOO 

NOW  USE  THESE  SS  TO  APPROXIMATE  THE  NEXT  BETAS 

CALL  CISPLN(BETA,S-E,THT,KJMX-2, 128,1) 

CALCULATING  THE  ARG ( DOMEGA/DZETA ) 

PB2=H?BN*FL0AT(N) 

DO  14  I=1,N2 

1  4  THT  ( I  )  =Ti-(T  (  I  )  -PHI  (  I  )  -PB2 

NON  FIND  THE  TIJ  AND  DRJ  COEFFICIENTS  CORRESPONDING 
ARG (DOMEGA/DZETA)  DATA 


DO  510  JP=1,M 

310  ZA<  JP)  =CMPLX(THT(2^<-JP-1  )  -THT(2*JP)  ) 

DO  307  NP=1,N 

307  ZA(NP) =CONJG(ZA(NP) ) 

.L  FFT(ZA-G  -  IWK  ) 

DU  308  JP=1,N 

308  ZA( JP) =CaNJG<ZA( JP) ) /ZNN 
CfV  L  SHUFL(N,ZA,ZCC) 

Zt>  :NP1  )  =ZA(  1  ) 

DO  311  NP=1,NB2P1 

2AJ (NP)=(C0NjG(ZA(NP2-NP) )+ZA(NP) )/2. 

3 . 1  ZA2(NP)  =ZHKC0N  JG(ZA(NP2-NP)  )-ZA(NP)  )  /2. 
Zi«  =ZN2N 

DO  312  NP=i,NS2Pl 
ZM=2N/ZH2N 

312  ZCC(NP)  =  (ZAl  (NP;+ZA2(NP)»ZUl)  /  2. 

Z’.-  -Zl'/ZW 

r  513  I=2,NB2 

ZH=ZW*ZW2N 

NP=NS2+I 

ZSB= ( ZAl ( NP2-NP ) f ZA2 ( NP2-NP ) *ZW) /2 . 

313  ZCC(NP)=CONJG(ZBB) 

I-.  -  (  i  )  =REAL  (  zee  (  1  )  ) 

ZCC(NP1 ) = .5*C0NJG(2A1 ( 1 )-2A2( 1 ) ) 

Di  (NPl )=REAL(ZCC(NP1 ) ) 

DC  314  NP=Z,N 
D  (NP)=2.^REAL(ZCC(NP) ) 

314  DR(NP)  =-2,-»AIMAG(2CC(NP)  ) 
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NOW  USE  DIJ  AND  DRJ  TO  GET  THE  NEXT  APPROXIMATION  OF 
ABS(D0MEGA/D2ETA) 

D0DZT0=AL0G(E<2) /PBN) 

Zi  Cd  )=CMPLX(0.  fO.  ) 

DfJ  302  JP  =  2,N 

302  ZCC( JP)=CMPLX(DR( JP)/2. .DK JP)/2. ) 

ZCC(NP1 ) =CMPLX(0. ,0. ) 

2-i  =  Zl/ZW2N 
DO  303  NP=i fNBZPI 

2m (NP) =2CC(NP)+C0NJG<2CC(NP2-NP) ) 

Zf>==ZW»ZW2N 

303  ZA2(NP)=2W*(ZCC(NP)-C0NJG(2CC(NP2-NP) ) > 

DO  304  NP=1 ,NB2P1 

304  ZA(NP)=ZA1 (NP)+ZI*ZA2(NP) 

DO  309  NP=2fNB2 

5i'S  ZA(NP2-NP) =C0NJG(ZA1 (NP) )+ZI*C0NJG(ZA2(NP) ) 

CALL  FFT(2A,Gf IWK ) 

CALL  SHUFL(NfZA,ZCC) 

DK  ( 1 ) =DCDZTO-REAL ( ZA  (  1  )  ) 

DO  20  K=1 ,N 

F(2»K-1) =EXP(DR( 1 )+REAL(2A(K ) ) ) 

;  2*K  )  =EXP(  DR  <  1  )+AIMAG<ZA<K  )  )  ) 

2(.  CONTINUE 

F;N2>1)=F(1) 

Gl  to  999 

900  WRITE(S,203)  ITrOM-DSMX 

2  '3  :-aRMAT( //IOXf  'DI  AND  DR  ITERATIONS  CONOERGED  AT  IT=  dlGr 
«  '  USING  CM  =  ',FS.3,'  .  MAXIMUM  ARC  LENG'H  ERROR  ^ 

iPEi2.3) 

CA_L  CISPLN(BETA,SfEfTHT,KJMX.2, 128, 1 ) 

X-'.TE(S,220) 

230  FORMAT (//4X,  '  K  '  ,  1 1 X ,  '  S  ' , 1 9X ,  'BETA',15X,  'PHI  ',//) 

DO  B9  K=1,N2 

GO  WR I TECS, 221)  K , E ( K ) , THT ( K ) , PH I ( K  ) 

221  FORMAT ( 15, 1P3E20.G) 

RETURN 

END 
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SUBROUTINE  SIMPSON ( ZO , ZN , A , B , M , N , ZI ) 
IMPLICIT  REAL(A-H,0-Y) , COMPLEX (Z) 

D' MENS  I ON  A(B5) ,3(G5) 

IP  (ZO.NE.ZN)  GO  TO  11 
ZI=0. 

RETURN 

CONTINUE 

ZM  =  CMPLX(2.'»-FL0AT(M)  ,0,  ) 

2k=(ZN-Z0)/ZM 

NM=N-1 

ZS20=CMPLX(A(N) ,B(N) ) 

ZSZN=CMPLX(A(N) ,B(N) ) 
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DO  10  J-l.NM 

ZSZO  =  ZSZO*ZO+CMPLX ( A ( N- J  > , B ( N- J ) > 
10  ZSZN  =  ZSZN*ZN+CMPLX(A(!M-J)  ,B(N-J)  ) 
Z  £0  =  CEXP(Z5Z0)+CEXP(ZSZN) 

Z i 1=CMPLX(0. ,0. ) 

Z;2=C«PLX(0. ,0. ) 

DC!  20  1  =  1,  MM 

ZZ=ZO+CMPLX( FLOAT < I ) ,0. )*2H 
ZSZZ=CMPLX( A(N) ,B<N) ) 

DO  30  J=1,NM 

30  ZSZZ=ZSZZ*ZZ+CMPLX(A(N-J) ,B<N-J) > 

(MOD(  1 ,2)  .NE.O)  GO  TO  31 
Z:2=ZI2+CEXP(2SZZ) 

GC)  TO  32 

31  ZI 1=ZI1+CEXP(2S2Z) 

32  CONTINUE 
20  CONTINUE 

Z: =ZH» <ZI0+2. *212+4. *ZI 1 )/3. 

R2TURN 

END 


SUBRCUT  :  NE  OMDZETA  (  20MST ,  ZTSTRT ,  A ,  B  ,  2MGA  ,  ZTAG5  ,  ZTANSR ,  N  ,  EPS  ,  i1 ,  NS  ) 
C 

C  THIS  SUBROUTINE  PERFORMS  THE  SAME  FUNCTION  FOR  ICIRC  =  1  AS  DOES 
C  O'^ETA  FCR  ICIRC  =  0,  NAMELY: 

C  -  c'CR  M  NOT  EQUAL  TO  1  OR  2,  IT  RETURNS  ZETA  (IN  ZTANSR)  FOR  A 
■:  GIUEN  UALUE  OF  SMALL  OMEGA  (IN  ZMGA),  USING  NEWTON-RAPHSON 

C  -  PCR  M  =  1,  IT  RETURNS  SMALL  OMEGA  (IN  ZMGA)  FOR  A  GIMEN  VALUE 
C  uF  ZE~A  (IN  ZTAGS) 

C  -  ‘=’CR  M  =  2,  IT  RETURNS  D(SMALL  OMEGA ) /D ( ZETA )  (IN  ZTANSR)  FOR 
C  GIVEN  VALUES  OF  ZETA  (IN  ZTAGS). 

C  THE  INTEGRATIONS  OF  D(SMALL  OMEGA ) /D ( ZETA )  ARE  DONE  IN  SUBROUTINE 
C  SIMPSON. 

C  ZOMST  IS  THE  VALUE  OF  SMALL  OMEGA  AT  THE  START  OF  THE  INTEGRATION. 

C  ZTSTRT  IS  THE  VALUE  OF  ZETA  AT  THE  START  OF  THE  INTEGRATION. 

C  NS  IS  THE  NUMBER  OF  STEPS  TO  BE  USED  BY  SUBROUTINE  SIMPSON, 
r.  N  IS  THE  NUMBER  OF  A  AND  B  COEFFICIENTS. 

IMPLICIT  REAL(A-H,0-Y) ,COMPLEX(Z) 

DIMENSION  A(S5) ,B(G5) 

Z-  =CMPL)<(  1  .  ,0.  ) 

2>--TA  =  ZTAGS 
)\M  =  N-1 
IT=1 

5  CONTINUE 

ZSMA=CMPLX( A(N) ,B(N) ) 

It  (M.EQ. 1 )  GO  TO  22 
IP  (M.EQ. 2)  GO  TO  40 

CALL  SIMPSON(ZTSTRT,ZETA,A,B,NS,N,ZSM) 

ZG=ZSM+ZOMST-ZMGA 
DO  10  J=1,NM 

10  2SMA=ZSMA*2ETA+CMPLX( A(N-J ) ,B(N-J) ) 

ZnG=CEXP(ZSMA) 


i:£TCLD  =  ZET  A 
ZETA=ZET0LD-ZG/2DG 

IF-  (CABS(ZETA-ZETOLD)  .LT.EPS)  GO  TO  20 

1  I  =IT+1 

IF  (CABS(ZETA) .GT. 1 . )  ZETA= . S*ZETA/CABS ( ZETA ) 
IF  (IT.LE.lOO)  GO  70  5 
i'i-  3 

kl-TURN 

20  ZTANSR=ZETA 
RETURN 

22  CALL  SIMPSONCZTSTRTtZETA^AfB.NS.N.ZMGA) 

Z.»1GA  =  2f1GA+20MST 
RETURN 

40  DO  41  J=1,NM 

41  ZSMA  =  2SMA*2ETA+CMPL)<(A(N-J  )  ,B(N-J)  ) 

2  I ANSR  =  CEXP(2SMA) 

RETURN 

END 


SUBROUT I NE  Z I SPLN ( Y , X , E , F , NP , IRTN , NRTN . NPD ) 


THIS  IS  A  COMPLEX-OALUED  UERSIDN  OF  CISPLN 

THIS  SUBROUTINE  FITS  A  CUBIC  SPLINE  To'^A  FUNCTION  Y{'''  =>-• 

NP  PAIRS  OF  POINTS.  THE  SECOND  DERI  VAT  I OF  THE  ^UNCTION 

=*PI.  .ir  NFB  =  I,  THE 

ASroED^NA^E^^'™°  "  Ordinate,  abscissa, abscissa. 
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IMPLICIT  C0MPLEX(A-H,0-Z) 

REAL  S I ZE  r  DF ,  DYDX  r  DA ,  DB .  DC  r  DL 

n^wlUlcTn!!  't(  ISO)  ,F(  ISO)  ,BDA(  ISO)  t  EM  (  1  SO  )  ,  H  (  1  Bn  ) 

DIMENSION  S( ISO) ,T( ISO) ,V( ISO) ,D( ISO) 

TPI  =  CMPLX(8.*ATAN(1 . ) ,0. ) 

THIS  SECTION  (ENTERED  WHEN  IRTN  =  1)  USES  THE  NP  PAIRS  OF  INPUT 
COORDINATES  X  AND  Y  TO  FIND  THE  COEFFICIENTS  OF  THE  SPLINE  FI~. 

~  CALLED  EM(KJ)  -  ARE  THE  SECOND  DERIVATIVES 

OF  THE  FUNCTiON, 


NPM  =  NP  -  1 
N  -  NP  +  1 

IF( IRTN.EG.Z)  GO  TO  20 
IF ( IRTN. EG. 3)  GO  TO  30 
IF ( IRTN. EG. 4)  GO  TO  40 
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DO  i  KJ  =  2rNP 

1  H(KJ)  =  X(KJ)  -  X(KJ-l) 

H<N)  =  H(2) 

DO  2  KJ  =  2rNP 

2  BDA(KJ)  =  K(KJ+1)/(H(KJ)+H(KJ+1)) 

Ell)  =  CMPLX(0.rO.) 

F( 1 )  =  CMPLXiO. ,0. ) 

Sa  )  =  CMPLX(  1  .  ,0.  ) 

DO  3  KJ  =  2-NPM 

DM  =  2.0  +  <1.0  -  BDA(KJ) )*E(KJ-1 ) 

E(KJ)  =  -BDA<KJ)/DN 
D  I  .<  J  )  = 

ft  S.0*( (Y(KJ+1 )-Y(KJ) )/H(KJ+l ) 
ft  -(Y(KJ)-Y(KJ-l) )/H(KJ) )/(H<KJ)+H(KJ+l) > 

5<KJ)  =  -< 1 .O-BDA(KJ) )ftS(K J-1 ) /DN 

3  F(KJ)  =  (D(K J )-< 1 .0-BDA(K J ) )*F{K J-1 ) ) /DN 
Y(N)  =  Y(2) 

IF(NPD.EQ.2)  Y<N)  =  Y<N)  +  TPI 
D(NP)  =  G.0ft(<Y(N)-Y(NP))/H<2) 
ft  -(Y<,MP)-Y<NPM)  )/H(NP)  )/(H(NP)+H(2)  ) 
iN’MM  =  NP  -  2 
T(NP)  =  CMPLX( 1 . ,0. ) 

0(NP,'  =  C^tPLXCO.  ,0.  ; 

DO  G  I  =  l.NPMM 
KJ  =  NP  -  I 

TUJ)  =  E<KJ  )*T(KJ+1 )  ft  S<:4J) 

0(KJ)  =  E(KJ)ftO(KJ+l )  +  F<KJ) 

G  CONTINUE 

F.'.  <NP)  =  (D(NP)-BDA(NP)ftO<2)-(l  .O-BDA(NP)  )*V(N?!'I/  )/ 
ft  (BDA(NP)*T(2)+( 1 .O-BDA(NP) )*T(NPM)  +2.0) 

DO  4  I  =  l.NPM 
K..  =  NP  -  I 

4  Ei^(KJ)  =  E(KJ  )ftEM(KJ  +  l  )  +  F(KJ)  +  S(1<J>*EM(NP) 

SIZE  =  AeS(X(NP)-X(l ) ) /lO.O 

RETURN 

C 

C  THIS  SECTION  (ENTERED  WHEN  IRTN  =  2)  RETURNS  NRTN  INTERPOLA 
r.  VALUES  OF  THE  ORDINATE  F  AT  ASSIGNED  VALUES  OF  THE  ABSCISSA 
C 

20  KJ  =  2 

DO  21  J  =  IrNRTN 
A  =  E(J) 

24  DA  =  REAL(A)  -  REAL(X<KJ)) 

D-  =  AIMAG(A)  -  AIi*!AG(X(K  J )  ) 

DA  =  DAftDA  +  DLftDL 

ns  =  REAL(A)  -  REAL (X(K J-1  )  ) 

DL  =  AII-AG(A>  -  AIMAG(X(KJ-1)  ) 

D^  =  DB*DB  +  DLftDL 

Dl.  =  REAL(X(KJ))  -  REAL(X(K  J-1  )  ) 

DL  =  Alf'AGCXCKJ)  )  -  AIMAG  (  X  <  K  J-1  )  ) 

DC  =  DC*DC  +  DLftDL 

IF  !  DA.l-E.DC.AND.DS.LE.DC)  go  TO  23 

KJ  =  KJ  +  1 

IFCKJ.LE.NP  )  GO  TO  24 
DF  =  ABS( A  -  X<NP) ) 
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l!^(DF.GT.SIZE)  W(?ITE(G,200)  J, E (  J )  ,NP, X(NP ) 

200  FORMAK //lOX,  'WARNING  -  ENTRY  IN  CISPLN  EXCEEDS  END  OF  BASE 

*  'ARRAY' ,/5X, 'E(  '  f  I3f  ' )  =  ',lP2EiB.8,'  EXCEEDS  X (',  15 =  ' 

*  ZEIG.S) 

D-  =  A8S(A  -  X( 1 ) ) 

IF<DF.LT. (-SIZE) )  WR I TE ( S , 20 1 ) J , E ( J ) , X ( 1 ) 

201  FORMAT(//iOX,  'WARNING  -  ENTRY  IN  CISPLN  IS  LESS  THAN  IRST' 

»  '  BASE  POINT', 

*  /5X,  'E(  ' , 13,  '  )  =  ',lP2ElG.a,'  IS  LESS  THAN  X(l)  =  15.8) 

KJ  =  NP 

23  DXA  =  X(KJ)  -  A 
DMB  =  A  -  X(KJ-l) 

CBA  =  DXA»DXA-*DXA 
CBS  =  DXB*DXB*DXB 

F(J)  =  <EM(KJ-1 )*CBA+EM(KJ)*CBB)/S.O/H(KJ) 

«  +  DXA*(Y(KJ-l)-EM(KJ-l)*H(KJ)*H(KJ)/G.O)/H(i<J) 

*  +  DXB*(Y(KJ)-EM(KJ)*H(KJ)*H(KJ)/G.O)/H(KJ> 

2'  CONTINUE 

RETURN 

THIS  SECTION  (ENTERED  WHEN  IRTN=3)  CONVERTS  THE  NP  PAIRS 
OF  INPUT  COORDINATES  X  AND  Y  TO  INTRINSIC  COORDINATES 
BETA  AND  S,  WHICH  ARE  STORED  AS  I 
BETA  *  REAL(E>,  S  *  REAL(F) 

NOTE  THAT  THE  SECOND  DERIVATIVES  OF  THIS  SPLINE  FIT  WERE 
ESTABLISHED  IN  A  PREVIOUS  CALL,  IN  WHICH  THE  FUNCTION  SMALL  OMEivi 
WAS  PASSED  IN  Y,  AND  THE  APPROXIMATE  ARCLENGTH  IN  K 

30  D(  1 )  =-H<2)*(EM(  1  )+EM(2>/2.  ) /3 .  ■►  ( Y  ( 2 ) -Y  ( 1  )  )  /H  ( 2 ) 

Di  31  KJ  =  2,NP 

31  D(KJ)=H(KJ)*<EM(KJ-1 )/2.+EM<KJ) ) /3 . + ( Y ( K J ) -Y ( K J- 1 ) ) /H ( K J ) 

DC  32  KJ=1 ,NP 

DL  =  ATAiM2(AIMAG(D(KJ)  )  ,REAL(D(KJ)  )  ) 

E<KJ)  =  CMPLX(DL,0. ) 

32  IF  (REAL(E(KJ) ) .LT.O. )  E(KJ)  =  E(KJ)+TPI 
DO  322  KK=2,NP 

IF( (REAL<E<KK) )-REAL<E(KK-l > ) ) .LT. ( -3 . 1 41 592G5 ) )E(KK ) =E(KK )+TPI 
322  CONTINUE 

E<NP)  =  E( 1 )  +  TPI 
F( 1 )  =  CMOLXCO. ,0. ) 

DYDX  =  AIMAG(D(  1  )  )/REAL(D(  n 

D-OTK  =  S(3RT(  1  .+DYDX*DYDX)*ABS<REAL(D(  1)  )  ) 

DC  33  KJ=2,NP 
DSDTO=DSDTK 

DVDX  =  AIMAG(D(KJ) )/REAL(D(KJ)  ) 


DYDX  =  AIMAG(D(KJ) )/REAL(D(KJ>  ) 

DSDTK  =  SQRTCl.  +  DYDX*DYDX ) ♦ABS ( REAL (  D ( !<  J  )  )  ) 
DL  =  ABS(X(KJ)-X(KJ-i ) ) 

33  F(KJ)  =  F(KJ-l)  +  DL*(DSDT0+DSDTK)/2. 

RC  TURN 


THIS  SECTION  (ENTERED  wHEN  IRTN  =  4)  CONCERTS  THE  NP  PAIRS 
OF  INPUT  COORDINATES  X  AND  Y  TO  INTRINSIC  COORDINATES 
SE~A  AND  5,  WHICH  ARE  STORED  AS! 

BETA  =  REAL(E)r  S  =  REAL(F) 

40  CONTINUE 

DO  41  KJ  =  2»NP 

41  H(KJ>  =  X(Kj;'  -  X(KJ-l) 

H  (  N )  =  H  (  2  ) 

DO  42  KJ  =  2.NP 

42  BDA(KJ)  =  H(K J+1 ) / (H(K J )+H(K J+1 ) ) 

E(  1  )  =  CMPLX(0. ,0. ) 

F(1 j  =  CMPLX(0. fO. ) 

S(  1  )  =  CMPLXC 1 . fO. ) 

DC!  43  KJ  =  2fNPM 

D\  =  2.  +  BDA(KJ)*E(KJ-1 ) 

E(KJ)  =  (BDA(KJ)  -  CMPLX( 1 . fO. ) )/DN 
D;KJ)  =  3.»(BDA(KJ:  XYCKJ  >-Y(K  J-1  )  > /H(K- / 

+  (C!^PLX(  ;  .  fO.  )  -  3DA<K  J  )  )♦»  V(K  J-f  1  ) -Y'K„  :  ' /H  '  K  J+ .  •  ! 

3(KJ)  =  -  3DA(KJ )*S(KJ-1 ) /DN 

43  !^(KJ)  =  (D(KJ)  -  BDA(KJ>*F(KJ-1  )  )/DN 
Y  (  N  )  =  Y  (  2  ) 

Ir(NPD.£Q.2)  Y(N)  =  Y(N)  *  TP  I 

D(NP)  =  3.*(BDA(NP)»(Y(NP)-Y<NP!>1)  ) /H^\P) 

*  (CNPLXd.  fO.  )-SDA(NP)  )»(Y(2)-Y(NP>)/H(2)  ) 

NPMM  =  MP  -  2 

T(NP)  =  CMPLX( 1 . fO. ) 

0(NP)  =  CMPLX(0.f0.) 

DO  4G  I  =  IfNPMM 
KJ  =  NP  -  I 

T(KJ)  =  E(KJ)*T<KJ+1 )  +  S(KJ) 

'.'(KJ)  =  £(KJ  )*'.'(K  J  +  1  )  +  F(KJ) 

43  CONTINUE 

Ei^(NP)  =  (D(NP)  +  <BDA(NP)-CMPLX(  1 .  fO.  )  )*'.'(2>-BDA(NP) 

«•  *'J(NPM)  )  /  (  (C«PLX(  1  .  fO.  )-BDA(NP)  )  *T  (  2  > +BDA  (  NP  >  ♦T  (  NP!*; ) 

♦  +  CMPLX(2. fO. ) ) 

DO  44  I  =  IfNPM 

KJ  =  NP  -  I 

44  EM(KJ)  =  E(K  J  )'»EM(K  J  +  1  )  +  F<KJ)  +  S(KJ)+EM(MP) 

SIZE  =  ABS(X(NP)-X( 1 ) ) /lO. 

D,i  47  KJ  =  IfNP 

47  D(KJ)  =  EM(KJ) 

DO  48  KJ  =  IfNP 

Du  =  ATAN2(AIMAG<D(KJ) ) fREAL(D(KJ) ) ) 

E;KJ)  =  CMPLX(DLfO.) 

48  IF(REAL(E(KJ) ) ,LT.O. )  E(KJ)  =  E<KJ)  +  TPI 
DO  422  KK  =  2fNP 

IF(  <REAL(E(KK  )  ) -REAL  (  E  (  KK-1  )  )  )  .  LT . -3 . 1 4 1 5S2B5  )  E  (  KK  >  =E(;<K  )■^^PI 
i22  CONTINUE 

E(NP)  =  Ed)  +  TPI 
Fd)  =  CMPLX(0.f0.) 


A  .sy- 


.s 
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D'DK  =  AIMAG(D(1) )/R£AL(D(l) ) 

DSDTK  =  SQRT( 1 .+DYDX*DYDK>*ABS(REAL(D( 1 ) ) ) 
DO  433  KJ  =  2,NP 
DSDTO  =  DSDTK 

DYDX  =  AIMAG(D(KJ) )/REAL(D(KJ) ) 

DSDTK  =  SQRT( 1 .+DYDX*DYDX)*ABS(REAL(D(KJ) ) ) 
D‘  =  ABS(X(KJ)-X(KJ-1 ) ) 

^33  F(KJ>  =  F(KJ-l)  +  DL*(DSDT0+DSDTK>/2. 

RETURN 

END 


SUBROUTINE  OMZTDR ( DR , D I , ZETAB , ZTA , 20MA , KMX r LMX , lOE) 

IMPLICIT  REAL( A-H,0-Y) ,COMPLEX(Z) 

THIS  SUBROUTINE  INTEGRATES  D(SMALL  OMEGA ) /D ( ZETA )  BY  THE 
TRAPEZOIDAL  RULE  TO  FIND  THE  GRID  IN  THE  SMALL-OMEGA  PLANE. 

IT  ,-,AKES  SOME  SMALL  ADJUSTMENTS  IN  THE  VALUES  OF  THE 
DERIVATIVES,  SO  AS  TO  MAKE  THE  RESULTS  PERIODIC. 

DIMENSION  ZTA( 10,40) ,ZDV( 10,40) ,ZSTR( 10) ,ZQMA( 10,40) 

DMENSION  DR(G5) ,DI (G5) 

Zr  =  CMPLX( 1 . ,0. ) 

DO  10  K  =  1,KMX 
DO  20  L  =  1,LMX 
ZTAGS  =  ZTA(L,K) 

CALL  0MD2ETA(2: ,21 ,DR,DI ,Z1 ,2TAGS,ZTANSR,G5, 1 .,2,1) 

2oV<L,K)  =  ZTANSR 
20  CONTINUE 
10  CONTINUE 

KMXH  =  <KMX+l)/2 
KM  =  KMXH  -  1 
KMXP  =  KMX  +  1 
KP  =  KMXH  +  1 

INTEGRATE  FIRST  ALONG  THE  PERIODIC  BOUNDARY: 

KASE  =  1 

27  IF( lOE.EQ. 1 )  GO  TO  11 

IF  lOE  EQUALS  0/1,  SMALL  OMEGA  =  -1.  IS  NOT/15  A  GRID  POINT 

CALL  OMDZETAIZI ,Z1 ,DR,DI ,21 ,ZETA8,ZTANSR,G5, 1 .,2,1) 

;KABE.EQ.2)  ZTANSR  =  ZR*ZTANSR 

ZOMA('_MX,KMXH)  =  -Z 1 +0 . 5-»  (  ZTANSR+ZDV  (  LMX  ,  KMXH  )  )  ^  (  ZTA  (  LMX  ,  KMXi-  ^ 
*  -  ZETAB) 

G''  TO  12 

i:  ZOMA(LMX,KMXH)  =  -Z1 
12  CONTINUE 

DO  25  J  =  1 ,KM 
K  *  KMXH  -  J 


.  •  o  -  ^ 
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Appendix  C 

DICTIONARY  OF  VARIABLES 


FORTRAN  SYMBOL 

ALGEBRAIC 

EQUIVALENT 

DEFINITION,  USE,  COMMENTS 

A(J).  B(J) 

"j  ■ 

Eq.  2-20,  2-21 

AK 

Eq.  2-23 

AKP 

Eq.  2-30 

AKQ 

- 

ALM 

A 

Eq.  2-31 

CAPK 

K  C^) 

Eq.  2-34 

CAPKPM 

kV^') 

Eq.  2-34 

DLT 

6 

A. 

Eq.  2-28 

DXIM,  DXIR 

Eq.  2-38,  4-4 

EX 

,/a-  f) 

^  r 

See  Figure  2 

EXINV 

^  ~  n 

* 

G,  H 

Gr,  H 

See  Figure  1 

HCPKPM 

j  K'(i') 

- 

OM 

— 

Relaxation  factor  used  in 
(p  ,  0  mapping 

FORTRAN  SYMBOL 


ALGEBRAIC 

EQUIVALENT 


DEFINITION,  USE,  COMMENTS 


ZPLUS,  ZMINUS 


Eq.  1-6 


Appendix  D 

LISTING  OF  METRIC  GENERATOR  PROGRAM 


LEVEL  21.7  (  DEC  72  ) 


OS/360  FORTRAN  H 


ISN  0002 
ISN  0003 
ISN  OOOA 
isri  0005 
ISi^  0006 
ISN  0007 


ISN  0008 
ISN  0009 
ISN  0010 
IS.I  0011 
ISN  0012 
ISN  0013 


ISN  0014 
ISM  0015 


'ILER  OPTIONS  -  NAME-  MAI N .OPT-02 . L I NECNT-60 , S IZE-OOOOK , 

SOURCE, EBCDIC. NOLI  ST. NODECK. LOAD, MAP. NOEO IT, ID. XREF 
C 

C  PROGRAM  CIMTVL  -  A  PROGRAM  TO  CALCULATE  THE  METRICS  OF  A  COORDINATE 
C  TRANSFORMATION  FOR  T'JRBOMACHINERV  CASCADES.  FOR  DOCUMENTATION . SEE : 

C 

C  0.  P.  NENNI  AND  W.  J.  RAE ,  EXPERIENCE  WITH  THE  DEVELOPMENT  OF 

C  AN  EULER  CODE  FOR  ROTOR  ROWS,  ASME  PAPER  83-GT-36.  MARCH  1983 

C  W.  0.  RAE,  A  COMPUTER  PROGRAM  FOR  THE  IVES  TRANSFORMATION  IN 

C  TURBOMACHINERY  CASCADES.  AFOSR  TR-81-0154.  ADA096416.  NOV  1980 

C  W.  J.  RAE,  MODIFICATIONS  OF  THE  I VES-L IUTERM02A  CONFORMAL- 

C  MAPPING  PROCEDURE  FOR  TURBOMACHINERY  CASCADES,  ASME  PAPER 

C  83-GT-116.  MARCH  1983 

C  W.J.  RAE.  REVISED  COMPUTER  PROGRAM  FOR  EVALUATING  THE  IVES 

C  TRANSFORMATION  IN  TURBOMACHINERY  CASCADES,  CALSPAN  CORPORATION 

C  REPORT  7177-A-l,  JULY  1983 

C 

C  THE  PROGRAM  READS  A  CARD  DECK  CONTAINING  THE  COORDINATES  X(K.L)  . 

C  Y(K,L)  ;  K  ■  l.KMX  i  i.  •  l.LMX,  AND  FINDS  BY  FINITE  DIFFERENCES 
C  THE  METRICS  OF  A  TRANSFORMATION  TO  A  RECTANGLE ( KS I , ETA > .  IN  THIS 
C  RECTANGLE.  THE  IllAGE  OF  THE  BLADE  SURFACE  LIES  ALONG  ONE  SIDE.  AND 
C  THE  IMAGES  OF  THE  POINTS  AT  INFINITY  LIE  AT  CORNERS  OF  THE  RECTANGLE 
C  (SEE  THE  REFERENCES  ABOVE).  THE  METRICS  ARE  WRITTEN  ON  TAPE  1. 

C 

implicit  REAL*8<A-H.O-2) 

DIMENSION  0(410,6) 

REAJ(5,100)  KMX.LMX 

100  F0R  'IAT(  2014  ) 

REA'.),  5, 101  )  SG 

101  FOR'UT(  8-10.4  ) 

C 

C  SG  IS  THE  SLANT  GAP  BETWEEN  BLADES 
C 

DO  10  K  •  l.KMX 
KMILMX  -  <  K-1 )«LMX 
LKA  •  KMILMX  1 
LKB  •  KMILMX  ♦  LMX 

10  READ(5.200)<<a(LK.5).Q(LK.6)),LK-LKA,LKB) 

200  FORMAT! 1P4E20. 13) 

C 

C  0(LK.5)  CONTAINS  XCK.L),  Q<LK.6>  CONTAINS  Y(K,L) 

C 

KM-.<MX-1 

LM-LMX-1 


—  -  SET  X  AND 


•KMX  AND 


ISN  0015 
ISn  0017 
ISM  0018 


ISH  0019 
ISM  0020 


LK-KM»LMX+I 
Q(  L.K.S  )-Q(  1,5) 
Q(LK,6)-Q< 1,6) 

—  FIND  X  AND  Y 


■1  AND  L-LMX 


Q(L;:X.5)-2DO*0(LM.5)-Q(LMX-2.3) 
a(L,'1X.6  )-2D0*O(  LM,5)-Q<LMX-2.6) 


AND  Y  AT 


■KMX  AND 


V.VVVV/T 


ISN  0021 
ISN  0022 
ISN  0023 
ISN  002« 
ISH  0025 

ISN  0026 
ISN  0027 

ISN  0023 
ISN  0029 
ISN  0030 
ISN  0031 
ISN  0032 
ISN  0033 

ISN  0034 
ISN  0036 
ISN  0033 
ISN  0040 


ISN  0042 
ISN  0043 
ISN  0044 
ISN  0045 
ISN  0046 

ISN  0047 
ISN  0049 


ISN  0061 
ISN  0052 
ISN  0053 
ISN  0054 
ISN  0055 
ISN  0056 
ISN  0057 


ISN  0059 
ISN  0059 
ISN  0060 
ISN  0061 
ISN  0062 
ISN  0063 


ISN  0064 
ISN  0065 
ISN  0066 
ISN  0067 
ISN  0063 


LK-<MX*LMX 

Q(  L:<.5  >«2D0«a<  LK-1 .5  )-Q(  LK-2.5  > 

Q<  LX,6  )«200*Q<  LK-1 .6 )-Q<  LK-2. 6  > 
TDELXI-400/0FLOAT(KM> 

TDElZT-2DO/DFLOAT<LM> 

C 

DO  520  K-l.KMX 
KM1LMX»( K-1 )*LMX 
C 

DO  510  L-l.LMX 
LK-:<M1LMX  +  L 
LPl »LK+1 
LMl-LK-1 
KPl^LK+LMX 
KMl-LK-LMX 
C 

IF<  <.EQ. 1  )  GO  TO  SOI 
IF( K.EQ.KMX  )  GO  TO  501 
IF( L  -EQ. 1  )  GO  TO  504 
IFIU.EO.LMX)  GO  TO  505 
C 

C  —  K-2  TO  KM  AND  L-2  TO  UM  (CENTERED  DIFFERENCES) 

C 

OXDXI-(Q( KPl .5 )-Q(KMI ,5 ) )/TDELXI 
DYDXI«<Q( KPl ,6  )-Q( KMl ,6  )  >/TOELXI 
DXD2T-<Q<LP1 .5)-Q<LMl,5) >/TDELZT 
0YDZT«(Q:LP1 .6)-Q(LMl,6) )/TDELZT 
GO  TO  509 
C 

SOI  IF( L .EQ.  1  )  GO  TO  510 
IFIL.EQ.LMX)  GO  TO  510 
C 

C  -  K«1  OR  KMX  AND  L«2  TO  LM 

C 

LK2-LMX+L 
LKKM-( KM-1  )*LMX  +  L 
DX0XI-( Q<  LK2,5  )-a( LKKM.S  )  )/TDELXI 
DYD-d-COI LK2,6>-Q( LKKM.S)  )/TDELXI 
OXOZT-( Q<  LP 1 .5  )-Q( LMl ,5  )  )/TDELZT 
0Y0ZT-<Q( LPl .5  )-Q( LMI .6 ) )/TDELZT 
GO  TO  509 
C 

C  —  L«1  AND  K-2  TO  KM 
C 

504  LKl-KMlLMX+1 
LK2-LK1*! 

LK3-LK2*! 

0X0ZT-< -300*a<LKl .5 )*400*Q<  LK2 . 5  ) -OC  LK3 , 5  )  )/TDEL2T 
DYD:T-( -3D0»Q<  LKl ,6  )  +  400*Q( L:<2.6  ) -Q<  LK3 , 6  )  ) /TDELZT 
GO  TO  506 
C 

C  — -  L-LMX  AND  K-2  TO  KM 
C 

505  LKS-< KMX-K  )*LMX  +  LM 
LK3-KM1LMX+LM 

DXDZT-(Q{ LKS.S )-a(LKB,5 ) ) /TDELZT 
DQ-SG 

IF( K.GT.KMXH  )  DQ--3G 


ISN  0070 


C 


DYDZT-(Qt  LKS.6)*0Q-Q(LKB.6) )/TDELZT 


ISN 

0071 

506  0X0;<I«(a(KPI  .6  )-a<KNI  .5  >  )/TDELXI 

ISN 

0072 

DVDXI-(Q(XP1 ,6)-Q(KMl .6) )/TDELXI 

ISN 

0073 

C 

c 

509  R0M1«(DXDXI*0V0ZT-0X0ZT«0VDXI > 

STORE  DERIVATIVES  IN  THE  ORDER  0KSI/DX,DKSI/DV,DETA/DX.0ETA/0V 

ISN 

0074 

Q<LT. 1  )-0Y0ZT/R0Ml 

IS.‘I 

0073 

Q<LK.2)— OXOZT/ROMl 

ISN 

0076 

a<LX.3)— OYOXI/ROHl 

ISN 

0077 

Q(LK.4)-  OXOXI/RDMl 

ISN 

0073 

510  CONTINUE 

ISN 

0079 

c 

c 

520  CONTINUE 

STORE  DERIVATIVES  ON  TAPE  1 

ISN 

0080 

c 

c 

LKMX-LMX*KMX 

WRITE! 1)  KMX.LMX.LKMX.I (Q( I . J ) . I - 1 , LKMX  > , J- 1 , 4 ) 

ISN 

0087 

c 

STOP 

ISN 

0088 

END 
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